For Reference 


NOT TO BE TAKEN FROM THIS ROOM 





For Reference 





NOT TO BE TAKEN FROM THIS ROOM 





DBOROODOOOONO 
The University of Alberta 
Printing Department 
Edmonton, Alberta 
ga Eg a re 


» “we ae 


UNIVERSITY OF ALBERTA 
LIBRARY 


Regulations Regarding Theses and Dissertations 


Typescript copies of theses and dissertations for Master's and Doctor's 
degrees deposited in the University of Alberta Library, as the official Copy of 
the Faculty of Graduate Studies, may be consulted in the Reference Reading Room 
only. 


A second copy is on deposit in the Department under whose supervision the 
work was done. Some Departments are willing to loan their copy to libraries, 
through the inter-library loan service of the University of Alberta Library. 


These theses and dissertations are to be used only with due regard to the 
rights of the author. Written permission of the author and of the Department 
must be obtained through the University of Alberta Library when extended passages 
are copied. When permission has been granted, acknowledgement must appear in the 
published work, 


This thesis or dissertation has been used in accordance with the above 
regulations by the persons listed below. The borrowing library is obligated to 
secure the signature of each user. 


Please sign below: 


Date Signature Institution 









imal cama 
Wyeiegupgnen eonoes.wensa 


a - wet oy = he 4 ‘ : 
- CHWYL r 
by /‘Bhwas. ; ye \ P 
noe ‘ ¥ 

a Tus eae a oe - 

a / 

a 
*, * 


teh am runaie a 
po OF cRapdare STUDIES 





THE UNIVERSITY OF ALBERTA 


AN ANALYSIS OF TRANSIENT GAS FLOW 


THROUGH POROUS MEDIA 


by 


EDWARD CHWYL 


A THESIS 
SUBMITTED TO THE FACULTY OF GRADUATE STUDIES 
IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE 
DEGREE OF MASTER OF SCIENCE IN 


PETROLEUM ENGINEERING 


FACULTY OF ENGINEERING 


DEPARTMENT OF CHEMICAL AND PETROLEUM ENGINEERING 


EDMONTON, ALBERTA 


FEBRUARY, 1968 


















ATAGSIA FO YTI@AAVING BHT. 


WOJT 2AD THSIGMAAT FO 2IeYJAVA UA 
ATGaM B2UOAOD HOUOAAT  . 


yd 
. IYWHD asAwaa 


 sismr a tet 

aatqura ZTAUGASD TO YTAUOAT. HT OT aay 
qHT dO% QTVSMAALUQGA BHT 10 TMSMITAIUG. 
Ul SOuaIDe YO ASTAAM IO saAaq 


UIASIUIONS MUAMOMTSG 


OUIAZAUIOWA FO vrauoas 
OMIATIMVIOUA MUSJOAT IS * ees 


UNIVERSITY OF ALBERTA 


FACULTY OF GRADUATE STUDIES 


The undersigned certify that they have read, and 
recommend to the Faculty of Graduate Studies for acceptance, 
a thesis entitled AN ANALYSIS OF TRANSIENT GAS FLOW THROUGH 
POROUS MEDIA submitted by Edward Chwyl in partial fulfilment 
of the requirements for the degree of Master of Science in 


Petroleum Engineering. 


igus | 








ADASEIA IO: YTLesavIuU 
esIaure aTAUGAD FO YTIUOAT a 


bus ,bsex sven yedt teddy ytisteo henpiexzebau: edT 
~2oasIqeoDs Ot eeibust2 etsubsx1D to ytivost eft oF 


HOUOHHT WOJY 24D TMURIGMAAT FO 2ISYIAMA WA belsisne eieed?- 
‘toromlitivt {ie6it1sq anit Lywdd Bbiswbha yd besttimdye AIGaM 2U0A0 


So » 







ak eonein2 to 1st32a6M toc esipsb edt. 10% esnemeztiups edt : 
v7 


ate ste te ates ry > -~ 


ABSTRACT 


It was initially thought that Darcy's law and the 
continuity equation Satisfactorily described transient iso- 
thermal gas flow through homogeneous and isotropic porous 
media. In the process of development, however, it was 
realized that two distinct phenomena, previously neglected, 
affected such flow. These were the inertial and the molecu- 
lar streaming effects. 

Although previous attempts have been made to pre- 
dict the effects of these phenomena on transient gas flow, 
no studies have been conducted where both effects have been 
considered simultaneously. In addition, solutions that 
account for either of these effects have generally been re- 
stricted to the case of a sealed external boundary and a 
constant mass producing rate. 

In order to verify the existence of inertial 
effects, the Forchheimer equation was derived for the linear 
horizontal system and extended to include plane radial flow. 

Existence of molecular streaming was justified by 
an extensive literature review on the subject. 

A mathematical model which accounted for both the 
inertial and the molecular streaming effects was developed 
for the linear horizontal and the plane radial porous systems. 


The Crank-Nicholson Implicit procedure, in conjunction with a 
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function iteration scheme, was employed to solve the linear 

case. It proved to be both stable and convergent for the 

range of variables considered. Numerical computation was 

carried out on an IBM 360/67 digital computer for the fol- 

lowing boundary conditions: 

1. Constant terminal rate with constant pressure at the 
external boundary. 

2. Constant terminal rate with a sealed external boundary. 

3. Constant terminal pressure with a sealed external 
boundary. 

4, Constant terminal pressure with constant pressure at 
the external boundary. 

The results obtained indicated that the inertial 
resistance coefficient affected transient gas flow more 
adversely than did molecular streaming. However, both 
Phenomena significantly affected the flow considered; there- 
fore, neither effect may be neglected without some prior 
consideration to determine whether or not it is in fact in- 


significant, 
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INTRODUCTION 


Although considerable effort has been directed to 
the study of transient isothermal gas flow through porous 
media, the condition of this topic remains incomplete. The 
major difficulty arises from nonlinearity of the partial dif- 
ferential equations that describe such flow. Consequently 
the only solutions available are approximate ones - graphical 
solutions, analogue solutions, approximate analytical solu- 
tions, and numerical solutions. 

Initially attempts were made to solve the non- 
linear partial differential equation resulting from a combina- 
tion of Darcy's law and the equation of continuity. In the 
process of development, however, it was realized that two 
distinct phenomena which were previously neglected signifi- 
cantly affected the physics of isothermal gas flow through 
porous media. These were the "inertial effects" and the 
"molecular streaming effects". 

Ever since the importance of these two effects 
was realized, there have been many attempts to predict the 
results each phenomenon would have on fluid flow through 
porous systems. Recent studies ‘t'2) have verified the fact 
that both these phenomena are simultaneously important in 
steady state laboratory core testing. 

It was the primary objective of this investigation 


to develop a workable mathematical model which would simul- 
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taneously incorporate these two phenomena and predict what 
effect they would have on transient isothermal gas flow 
through finite porous systems. 

A secondary objective of this thesis was to validate 
the theory underlying visco-inertial gas flow and, to a lesser 


(3) 


extent, molecular streaming. Hubbert succeeded in deriving 
Darcy's law from the Navier-Stokes equation of motion. It 

was felt that a similar procedure could be extended to theo- 
retically verify the validity of Forchheimer's equation. The 


importance of molecular streaming could be established by an 


extensive literature survey on the subject. 
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LITERATURE REVIEW 


A. Viscous and Visco-Inertial Flow Equations 


(4) (5) 


Hubbert quotes Henry Darcy to have experi- 
mentally established that the downward flow of water through 
sand filters was governed by the following equation: 

q) =dexix analy, (1) 

H 

This was subsequently called Darcy's law. Ever since this 
pioneer work, attempts have been made to theoretically derive 
Darcy's empirical equation and to employ variables with more 
general physical meaning. Numerous procedures had been fol- 
lowed, including capillary tube analogies, hydraulic radius 
theories, drag theories, turbulent flow analogies, and di- 
rect derivation from the Navier-Stokes equation of motion. 
In view of the extensive literature published on this subject, 
the following literature review has been limited to procedures 
involving derivations from the Navier-Stokes equation. 

Using the Navier-Stokes equation of motion as a 
basis, Hubbert ‘°) derived Darcy's Law which was valid for any 
macroscopically homogeneous and isotropic porous medium. 
Treating viscosity as a constant and neglecting inertial 
forces, he integrated each term of the Navier-Stokes equation 


over the volume element, AV, and thus justified neglecting 


the divergence term (v-v) for gas flow as well as liquid flow. 
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He concluded that Darcy's law describing gas flow was identi- 
cal to that for a liquid "provided that the flow behavior of 
a gas in small pore spaces, other than expansion, is similar 


to that of a liquid". His final form of Darcy's law was 


2 
Yo in rll 


q = ——__——_——_—— (g - — grad P ) (2) 
Hu 9) 


Warren A. Ha11'”) derived Darcy's law for an in- 
compressible fluid flowing through anisotropic porous media. 
Also considering viscosity as a constant and inertial forces 
as negligible, he selected a volume element of fluid in. te 
porous medium "which is small with respect to the gross di- 
mensions of the system, yet large enough that the surface 
area of the matrix can be considered uniformly distributed 
throughout the element". He summed the body forces (pressure 
and gravity) and the surface forces (viscous) acting on this 


fluid element and obtained his version of Darcy's law. 


1 
gies K- (VP + pVgz) (3) 


- 
ned 


where permeability, K, is considered to be a tensor, and 
(gz) is the potential of the gravitational force. When the 
medium is isotropic, permeability becomes a scalar quantity, 
in which case the dot product degenerates to a scalar pro- 
duct and Equation (3) may be rewritten as: 

k 


G = - — (VP + pVgz) (4) 
u 
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S. Irmay considered one-dimensional flow of an 
incompressible fluid through an homogeneous and isotropic 
medium. Employing statistical methods on the Navier-Stokes 
equation, he succeeded in obtaining an equation that. de- 
scribed flow around spherical particles: 
(3,) = aa, + ba,” + cd, (t) (5) 
5 x x x 


Through dimensional analogies the constants a, b, 


and c were estimated to be: 


2 
N,, (1-4 )v 

a = See aT 2 5 . (6a) 
og ( p) 
No. (L=6) 

eee $4 ____ (6b) 
> g(d_.) 
i 

wea (6c) 
og 


Irmay stated that Polubarinova-Kochina (1952) showed the 
term, cq, (t), to be insignificant with respect to the other 
two terms. Equation (5) thus reduced to an. equation consist- 


(9) 


ing of two terms, as was suggested by Forchheimer 


(J,) =) Legg + ba, (7) 
x 


If, in addition, the quadratic term is neglected, a Darcy~ 


type equation results. This is illustrated by Equation (8). 


igeneme lage (8) 
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Sa Whitaker (19) | using Cartesian coordinates and an 


integrating technique similar to that of Huppae eels, derived 
Darcy's law for the flow of an incompressible fluid through 
an anisotropic but homogeneous porous medium. Introducing 
Euclidean shifters he was able to extend his derivation to 
general curvilinear coordinates, but for isotropic porous 


systems only. The generally valid form of Darcy's law as 


stated by Whitaker is: 


k 
q=- — (grad P - 0g) (9) 
u 


ae credit Houpeurt (1959) with 


Rowan and Clegg 
deriving a Forchheimer-tvpe flow equation for an idealized 
porous system in which each channel was approximated by a 
series of truncated cones. 

Except for these works of Irmay and Houpeurt, 
very little attention has been directed to the theoretical 
derivation of Forchheimer's equation. The alternate approach 
was to experimentally justify the existence of the quadratic 


term that is present in visco-inertial flow. The linear 


Forchheimer equation for horizontal flow may be written as: 


oP 2 
——— —= Fouqdy + Fy 0 (q,) (10) 
ox iy 
where - = is a "macroscopic" pressure gradient in the 


direction of flow. 
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If the truncated form of Equation (10) is compared 
with the horizontal Darcy equation for linear flow, it be- 
comes obvious that the viscous resistance coefficient, For 
is equivalent to the inverse of the permeability. 


(13) conducted flow tests on a num- 


Green and Duwez 
ber of sintered cores and deduced that gas flow through these 
cores could be described by an equation idential to Equation 
610) PheBoth FE and Fo were constants describing the porous 
medium and independent of the flowing fluid. Using the 
ideal gas law, the authors integrated Equation (10) over the 


length of the core sample and obtained the following equation: 


2 





=AP 2RT ZRT 2 
See) wee) WW (11) 
a x b x 
L gM GoM 
port 
where Bin ae pie 
Ergun 14) illustrated the validity of an expression 


similar to Equation (10) by showing that a plot of 


eAP 1 


L evs. (oq,) yielded a straight line relationship. 


(45) ave 
Therefore, he concluded that an equation similar to Equation 
(10) served as a satisfactory correlation for his empirical 
data. 

After conducting flow tests on twenty-four samples 
of sandstones, limestones, and dolomites, Cornell and Katz (19) 


were able to conclude that the phenomenon of flow of the four 


gases that they used could best be described by an equation 
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identical to Equation (11). 


(16) conducted flow tests on 


Greenberg and Weger 
samples of sintered metals and confirmed the results of Cor- 
nell and Katz. In addition, they showed that the term Re 


was a strong function of temperature whereas the term Fy was 


in fact independent of temperature. 


B. Phenomenon of Molecular Streaming 


In data published by Fancher, Lewis, and Barnes 
in 1933, there was a consistent discrepancy between the per- 
meability of oil sands to water and the permeability of oil 
sands to air. Prompted by this discrepancy, Klinkenberg ‘44! 
investigated the validity of the assumption that the per- 
meability of a porous medium is independent of the nature of 
the fluid. In what is now considered the classical paper on 
molecular streaming of gases in petroleum reservoirs, Klinken- 
berg concluded that the permeability to a fluid was indeed a 
function of the nature of that fluid. He approximated the 
porous medium by a bundle of capillaries, one-third of which 
were oriented in each of the coordinate directions. By de- 
riving a modified Poiseuille's equation which accounted for 
gas slippage in this model of the porous medium, and compar- 
ing the equation obtained with the integrated Darcy equation, 
Klinkenberg was able to conclude that the relationship be- 


tween the apparent permeability to a gas and the absolute 
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(liquid) permeability was 


Ace a), 


Kk. = k(l + geet Steves (12) 


He reasonedthat because the mean free path, i is an in- 


ave’ 
verse function of the mean pressure, Equation (12) could be 


rewritten as 


k = k(l+—) CES) 


where b is a lumped constant and is a characteristic of the 
porous medium and the particular flowing gas. Klinkenberg then 
experimentally demonstrated that Equation (13) was applicable 


to porous media in general. 


(18) (13)) 


Calhoun and Yuster and Rose concluded that 
the apparent permeability of a given porous medium to different 
gases is equivalent if the temperature is such that the mean 
free paths of the different gas molecules is equivalent. This, 
in essence, is a restatement of Equation (12) which states that 
the apparent permeability to a gas is directly proportional to 
the mean free path of that gas. Rose derived an expression 


which predicted the mean free path length at any mean pressure 


and temperature: 


r = 2.13 — RT/M (14) 
ave 


cette were among the first to 


Collins and Crawfor 
incorporate the Klinkenberg slippage correction term into an 


equation describing transient gas flow. They modified 
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Darcy's law for horizontal gas flow, Equation (15), by sub- 
stituting the Klinkenberg correction term for the apparent 


permeability. The result was Equation (16). 


k 

q = .-i-2- grad P (15) 
u 

% k  P+b 

q = --— (—) grad P (16) 
u P 


When they combined Equation (16) with the continuity equation 
and the ideal gas law, they obtained the following partial 


differential equation: 


5 2 2ud 3 (P+b) 
Vv (P+b) = — —— is) 
k Mh of 


The authors then instituted a change of the dependent vari- 
ables by setting P* equal to (P+b), and obtained a partial 
differential equation identical in form to one in which the 
slippage CeereckiGn term was not accounted for. This equation 
was 2ud 9 (P*) 


Ge(*) ~~ Sop —3 (18) 
k dt 





They therefore concluded that the solution to Equation (18) 
would be the same as the solution to the Darcy-Continuity 
equation, provided that the dependent variable was redefined 
in the manner specified. 


(21) 


Stewart and Owens experimentally demonstrated 
the reduction of gas slippage with increasing inertial effects. 


Plug flow was characterized by the complete disappearance of 
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the slippage phenomenon. The apparent permeability for plug 
flow was therefore independent of the mean flowing pressure 


of the gas and dependent solely on the Reynold's number. 


C. Solutions to Transient Isothermal Gas Flow Equations 


Various forms of the partial differential equation 
describing transient gas flow through porous media have been 
solved by different approximate procedures. Discussion, how- 
ever, will be restricted to the semi-analytical methods and 
the numerical methods. The cases treated are for isothermal 
horizontal flow of gases exhibiting constant viscosity and 
gas compressibility, unless otherwise specified. The porous 
media considered are homogeneous and isotropic. 

The basic partial differential equation describing 


horizontal transient Darcy-type flow of ideal gases through 


porous media, as stated by Bruce et at oe) is 
2ou oP 
aaGe pee) eee (19) 
Kkysat 


For flow in the x-direction, Equation (19) becomes 





3p? 2ou oP 
5 Sai Somer es (20) 
ox k ot 


For plane radial flow, Equation (19) becomes 


3 2p? ir 3p? 2ou dP 


5 +- —- 2 = o——— — (21) 
ox aa gi ad Ke ace: 
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The authors numerically solved Equations (20) and (21) for 
the case of a constant production flow rate of an ideal gas 
from a reservoir with a sealed external boundary. They em- 
ployed both the Forward Difference Explicit and the Crank- 
Nicholson Implicit procedures. (For a discussion of these 
finite difference procedures, refer to Lapidus ‘2?) ,) The 
solutions, presented in graphical form for various values of 
dimensionless parameters, were compared with a laboratory 
study of gas depletion in a linear system having the same 
boundary conditions. 


Jenkins and Aronofsky ‘74? 


presented a numerical 
technique for solving a radial system identical to that of 
Bruce et al. They redefined the independent variable as 

— = ln r, and used the Forward Difference Explicit procedure 
to solve for pressures at all the grid points except near 
the wellbore. At the wellbore they employed a "steady state 
core approximation" in order to avoid excessive computer 
time. 


Aronofsky and Porter '29) 


extended the above pro- 
cedure to include a real gas with varying viscosity and gas 
compressibility. The partial differential equation that they 
obtained to describe such flow was the following: 

k 30 


V(o —~ VP) = $6 — (22) 
U at 
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where 
vols (c5P + C3) 
2 eae C4P) 
P.2, Pp 
(9) = —<<<=—=<= << 
P ba 
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The terms Cor C3, and Cc, are constants which characterize 
the gas flowing through the porous system. 
The equation that they solved was the expanded 
radial form of Equation (22), namely: 
9°P a~P aP 
[B) to—aelRgdP) = F,(P) 45 2475 
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(23) 


where F,(P), F,(P), and F, (P) are functions of the character- 
istic values of the real gas at the flowing pressure of the 
system. 

The first attempt to solve an equation describing 
transient gas flow corrected for molecular streaming was 


ick an He utilized Darcy's law, the equa- 


made by Aronofsky 
tion of continuity, the ideal gas law, and the Klinkenberg 
slippage correction term to obtain the following partial 
differential equation: 


ou oP 
7+] (p+p) vP| = —— (24) 
kinet 


He then solved Equation (24) for the constant terminal pressure 


case with a sealed external boundary by using the Forward 


- ££ + 


. ’ net Ww oe aeteiie: 

(0+ haan Pee ae 

(4,9 e “ ‘eg a. a 
. ) ty Atel temas 


e.| : 
_ Sai ae 
s ae 





















S 
&- 
4 


aa Hake? 
esixzssosxsdo Asie arneruaee ot6 49 bas Edi wets 


-meseye euoxog orit a HF . vies 

ef a PY 

bebnsqxs ont 26W bevioe ysds tees soiseupe 
tc wit 
: \lemea ($8) rag to a 


= S 

(€8) . 3898 we = see foe. 
Gee 263-0 
~retosisio edit to eisieiglanal e158 ygtinen DT yt 
eit to sxveesiq paiwolt ed? 36 me: thex wth uLs 
wins ia a a 

paidizoasb noitsups as svioe ot aha Aa enh | 
esw prtimsorse isivoslom xot bst29x102 ¥ _snaies 
~sups eft ,wsl e'yousd bestlizu ef “ma 2 


prednednila od bas ,wsl esp Leebi eats antes Mi 
Isitxeq eaiwosted exit atete Ae oex 109 9 it 


oe 
au..; 
2 





~, 


a? 


Difference Explicit procedure. This work was later examined 


experimentally by Wallick and Aronofsky ‘*’) , 


The experimental 
work showed excellent agreement with the theoretical pre- 
dictions except for the initial transient period where flow 
rates were the highest and inertial effects were the most 
Significant. 


Kidder ‘28? 


presented an analytical solution to the 
problem of transient flow of a gas through one-dimensional 
semi-infinite porous media with a constant pressure at the 
producing face. He employed a perturbation technique 


(carried out to include second-order terms) and obtained a 


solution to the following partial differential equation: 


Pe) dP gu dP 
2 WEG ies (25) 
ox ox k at 
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Al-Hussainy et anilee considered the effect of a 
pressure dependent viscosity and compressibility factor on 
Darcy-type flow of real gases through porous media. By 

eliminating density from Equation (22),the authors obtained 


the equation 
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The above equation was transformed into a form similar to the 


diffusivity equation, but in addition it relaxed the assump- 


tion that pressure gradients were small. The result was the 
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following diffusivity-type equation: 











5 du(P)c_ (P) om(P) 
Can (hye (27) 
k at 
where 
Pp P 
m(P) si 2 —————._ dP 
P, u (P) z (P) 
l do 
e{g) = ee 
0 AP 
In the radial system, Equation (27) becomes: 
3°m(P) 1 am(P) ¢u(P)c_(P) am(P) 
5 2 era = eee (28) 
or ror k dt 


The authors concluded that the solutions obtained 
to Equation (28) in terms of the pseudo-pressure drop, m(P), 
were analagous to the solutions that van Everdingen and 


Hurst (29) 


obtained to a similar equation, but with a dimen- 
sionless pressure drop as the dependent variable. 

Tek et We Get combined the radial Forchheimer 
equation, which was similar to Equation (29), with the radial 
continuity equation, and numerically solved the resulting 


partial differential equation using the Backward Difference 


Implicit procedure. 
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The solutions were restricted to the case of a constant 
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terminal production rate and a sealed external boundary. By 
varying the inertial coefficient the authors were able to 
study its effect on wellbore pressures, drainage radii, and 
isochronal testing. 

Swift and Kiew (32) numerically and semi-analytically 
solved the partial differential equation resulting froma 
combination of the radial Forchheimer equation with the 
radial continuity equation. By introducing a correction term, 
6(r,Q), to compensate for non-Darcy flow, they were able to 
obtain a partial differential equation very similar to one 


describing Darcy flow, namely: 
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Equation (30) was solved by a numerical technique 


eee 


Similar to the one employed by Bruce eta for a reser- 


voir with a constant producing rate and a sealed external 
boundary. The authors expressed the Forchheimer equation 


as: 
dP 
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Rowan and Clegg '34) 


semi-analytically solved Forch- 
heimer-type radial gas flow through bounded and infinite 
reservoirs producing at both constant rates and constant 
pressures. They employed a radius of disturbance which 
"moved outward from the wellbore and was such that beyond 

it the formation was undisturbed". Because a mean value for 
(9P/dt) within this zone of disturbance was used, the Forch- 
heimer and the continuity equations could be rewritten as 
ordinary differential equations and therefore integrated be- 
tween the appropriate limits. The semi-analytical solu- 
tions were compared with published results for quadratic 
flow obtained on digital and analogue computers. This pro- 
cedure was essentially an extension of an earlier work by 


(35) 


Rowan and Clegg where they considered Darcy-type flow in 


radial porous systems. 
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THEORY 


Isothermal laminar flow of a Newtonian fluid 
through any continuous flow channel can be described by the 
microscopic Navier-Stokes equation of motion, which is 
rigorously derived by Bird, Stewart, and Lightfoot (36) , 

The equation can be written in vector notation, ona per 


volume basis, as: 


t. T Dv U 2+ aR 
g-- grad Pp = — - ~|Vv+— V(Vev) S24 
0 Dt ie) 3 


Since complete geometrical characterization of a 
porous medium has not yet been accomplished, it is most dif- 
ficult to accurately describe flow through a porous system. - 
Geometrical parameters of a porous medium are based on volu- 
metric averages and, as a result, average flow values should 
be employed to describe the physics of flow through a porous 
System. Thus, an acceptable procedure would be to integrate 
the microscopic variables in the Navier-Stokes equation over 
some macroscopic volume element. Ideally, this element 
should be large enough so that the average quantities thus 
obtained generally describe the porous medium and the flow 
through it, and yet small enough so that the flow equations 
can still be considered continuous. Microscopic quantities 
would therefore be converted into measurable macroscopic 


quantities. These two frames of reference are best described 
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by Hubbert '37) , He states that the macroscopic scale "is 


large as compared with the grain or pore size of the porous 
solid and the flow of a fluid through a porous solid is seen 
as a continuous phenomenon in space". The microscopic scale, 
on the contrary, is "commensurate with the grain or pore size 
of the solid, but still large as compared with the molecular 
dimensions or motional irregularities due to Brownian or 
molecular movements." 

For the sake of brevity, the remaining theory has 
been confined to the horizontal linear and the plane radial 
geometries, which are illustrated by Figures 1 and 2 res- 
pectively. In particular the significance of inertial 


effects and molecular streaming effects are considered. 


A. Linear System 


For linear flow in the x-direction, Equation (32) 


can be written as: 
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Figure 2. Radial Geometry (A&é = 0.25) 
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If flow is further assumed to be horizontal and 


steady, Equation (33) becomes: 
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Since the Forchheimer equation is valid for steady 
state flow which accounts for viscous and inertial effects, 
it appears reasonable that it could be obtained from Equation 
(34), which describes similar flow but on the microscopic 
scale. 

This can in fact be accomplished, as indicated in 
Appendix A, by applying an appropriate integration scheme to 
Equation (34). This procedure yields the linear form of the 
Forchheimer equation, namely: 

dP 
- ei See ods te pa, |a,| (35) 
where (d3P/dx) is the macroscopic pressure gradient. 

If this equation is multiplied through by density 
and then density is eliminated from the left hand side by 


employing the equation of state, 


where Cc = 

















. v6 ve ‘ 
-- ee a v) = 

> —. Y 
ae ; 2 (9 oe 5 

| ¢ 
R >" 6 -Vv ¢ 2 
~ f “* 
a & cmebebe —* 4 f e pe + —— + 
; xt He € 5 ys 
c. * ¥ . 


wo 


{haste ot Dilay ef agktbype vent: an dioxc ona sone . 
,atoet?ts [sitren: bons essoelvY oe etapoo8 dokdw vot a z 
gadi © Idsnoesst euneqge 
: chal’ om as 


otqo0o7veu0lm sift an ted wold’ wehbe eedi zeae fo iv 


*. 


/f 
3) 
iG 
te 


*. 


yi 


i 


. 


wt mso et Lat Hy 


st 


bersotbai es ,bedeliqnoods o@ F983 


i 
TZ 
manos {9 isY¥uet: TE - CO toEdsS ns eaiytags veh A ute ut ry 
a : ah 
niot isecil aid ebleby exubedoxg eld? i(88) 4 otter 


2 . 
¢ 


: yiemsa ,aot _— 


Pu 4a sg, Lee al -~) 


G 7 ba 


i 
; 4 : ipee 
Y7feash ya NpwotAy Dbedig Jf um ex ddksaupe'e 





Pa coke a 


a more workable form results, viz., 


3p? 2 u | 
-— = — | — (pq) + F, (pq) |oq.| 87) 
a G k x b x x 
fe) a 
Wi 
where —_— = FE. 
k 
a 


Since the term kK. in this equation denotes the 
apparent permeability, which is a function of the flowing 
pressure and the geometry of the porous structure, it is 
appropriate to substitute for it by means of the following 
equation: 

b 
k = k(l + —) (38) 
a 
3 
When this substitution is made a modified Forchheimer equa- 


tion, which accounts for the molecular streaming of gases, 


3p? 2 u P 
= — = —| — (—)(leg,) + Fi leg.) lea, |] (39) 
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results: 


To describe transient gas flow, the above flow 
equation must be coupled with an appropriate equation of 
continuity, namely: 


2 
P) st dP 


ek (eq) aie, ee ee (40) 
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The equation of continuity is valid for all times 
whereas the modified Forchheimer equation, which describes 


steady state flow, is valid at any particular instant of 
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time. A combination of the two equations will therefore 
amount to a series of steady state approximations to the 
transient flow problem, Once Equation (39) and (40) are 
coupled with an appropriate set of boundary conditions, then 
linear, horizontal, isothermal, transient flow of gases is 
fully described. 

In order to make the approach general for a wide 
range of physical properties and flow conditions, the modi- 
fied Forchheimer equation and the equation of continuity may 


be rewritten in terms of dimensionless parameters as: 


ap? p O . 
-—— = 2(—) (5q.) + 4B($q_) |bal| (41) 
ax P+b % - a 
and 
abt 1k, ap 
SS iets Wa fr cae (42) 
ax 2Pevats 


If these two equations are combined according to a 
procedure outlined in Appendix B, a second-order nonlinear 
partial differential equation results which may be written 


in the following manner: 
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The dimensionless flow rate per unit area, eds 
appears in both these nonlinear terms and is related to the 
dimensionless pressure distribution by Equation (41). 

To completely describe transient horizontal linear 
gas flow through a porous system, Equations (43) and (41) 
must be coupled with an appropriate ens of boundary condi- 
tions. The four most common sets of boundary conditions for 


a finite linear porous system are summarized in Table l. 


B. Radial System 


The Forchheimer equation for plane radial flow may 
be expressed as: 
3p- 2 u 
oe vl oes ae 9 LT AD inate an i fone Ba (22° Baa (44) 


or fe. k 
Oo a 


This equation, which was developed in Appendix A, is identi- 
cal to the linear Forchheimer equation except that the in- 
dependent variable (x) is replaced by (r). 

If Klinkenberg's slippage correction term is intro- 
duced into Equation (44), a radial Forchheimer equation cor- 
rected for molecular streaming results. This equation can 


be written as: 


3p? 2 u P | | 
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Also needed to describe transient flow behavior is the radial 


equation of continuity, namely: 
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Introducing dimensionless parameters into Equations (45) and 


(46), they may be rewritten as 
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Equations (47) and (48) may be combined according to the 
procedure outlined in Appendix B, the result being the 


following second-order nonlinear partial differential equa- 


tion: 
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In order to permit closely spaced increments near 
the wellbore and still retain equally spaced increments of 
the space variable, the transformation, & = ln r, may be 


employed. Under this transformation Equation (49) becomes: 


C(u,u,)u,, + D(u,u,)uU, + E(u,u,) = Nee (50) 
where A oe pe 
P(P+b)e 78 
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The dimensionless flow rate per unit area, (oq, )y 
once again appears in all the nonlinear terms. It must be 
related to the dimensionless pressure distribution by the 
Forchheimer equation rewritten in terms of the new space 


variable. This equation may be expressed as: 


3P- 2P (pq,) Si otal ten OS 
~-—— = ———— + 48(0g,) |p, le (51) 
J & (P+b) 


To completely describe plane radial transient gas 
flow through a porous system, Equations (50) and (51) must 


be accompanied by an appropriate set of boundary conditions. 
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The four most common sets of boundary conditions for a finite 


radial porous system are summarized in Table 2. 
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SOLVING THE MATHEMATICAL MODEL 


If it is conceded that Equations (43) and (41) and 
Equations (50) and (51), along with an appropriate set of 
boundary conditions, satisfactorily describe the flow of gases 
in their respective porous systems, the task remains to obtain 
reasonable solutions to these formulated equations. The main 
objective of such solutions is to predict the effect that the 
dimensionless slip and inertial resistance coefficients ex- 


hibit on transient gas flow through finite porous media. 


A. The Procedure Employed 


Equations (43) and (50) can be classified as second- 
order, nonlinear, parabolic-type partial differential equa- 
tions. Primarily because of their nonlinearity, these equa- 
tions do not lend themselves to analytical solution and con- 
sequently must be solved by means of a numerical technique. 

The technique must be stable and it must converge to solu- 
tions that can be considered correct. Since it is necessary 
to couple Equations (43) and (50) with Equations (41) and (51), 
respectively, the complexity of the problem increases. 


Hmesor ce 


classifies finite difference approximation 
procedures as either explicit or implicit. Explicit proce- 
dures yield solutions directly in. terms of known quantities, 


but stability problems arise if the mesh size ratio is not 
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properly chosen. Implicit procedures involve increased com- 

putational complexity but overcome this stability problem. 

In order to minimize problems of stability, implicit proce- 

dures were favored in general and the Crank-Nicholson proce- 

dure was chosen in particular. Specifically, the following 
discrete approximations were employed in the numerical 
technique: 

1. A second-order correct central difference approximation 
on the space derivatives in the partial differential 
equation. 

2. A second-order correct central difference approximation 
on all Neumann boundary conditions, 

3. A third-order correct approximation on the space deriva- 
tive in the Forchheimer equation. 

These finite difference approximations yielded a 
system of nonlinear algebraic equations which may be written 


in matrix notation as 
f > i: => 
a ar cg Py eal Rac bin (52) 
The square root inversion technique, modified to accommodate 
nonsymmetric matrices, was employed to invert the matrix. The 


solution vector was then obtained according to Equation (53) 
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From the values of the pressures thus obtained, the flow rates 
at time level (m+l) were evaluated. As a result, the non- 
linear terms could also be calculated. 

In order to ensure that the flow rates and the 
pressure distributions agreed within a reasonable degree of 
accuracy, a material balance was performed for each boundary 
condition. 

The computational scheme employed can best be il- 
lustrated by means of a general flow diagram as shown in 
Figure 3. A detailed copy of the computer program for the 


second set of boundary conditions is presented in Appendix E. 


By Stability and Convergence 


As has been previously stated, the numerical pro- 
cedure must first of all be stable, and secondly, it must 
yield a solution that is "reasonably close" to the exact 
solution of this system of equations. 

A numerical procedure can be termed "stable" if an 
error introduced at the outset of an iteration scheme does 
not magnify itself as iteration continues. Instability prob- 
lems, which would have eventually swamped the true solution, 
were not encountered in the numerical technique employed. 

Normally a numerical solution can be termed "con- 
vergent" if it approaches the exact solution as the grid 
spacing approaches zero. Since the exact solution to the 


problem in question cannot be deduced, a different con- 
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Figure 3. General Flow Diagram 
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vergence criterion was necessary. Consequently convergence 
was assumed when the results obtained at one grid spacing 
agreed within a certain "allowable error" to those obtained 
at a grid spacing that was one-half the previous one. Both 
the dimensionless space increment (Ax) and the dimensionless 
time increment (At)) were checked. This was not an absolute 
check on convergence since the numerical results could still 
differ from the exact solution depending on the rate of 
convergence. However, for want of a better method, this 


convergence criterion was employed. 
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RESULTS 


The Forchheimer equation describing linear, hori- 
zontal, steady gas flow was derived from the Navier-Stokes 
equation of motion by employing a procedure similar to that 


of Robbertic |. 


The-result wassEquation (37). MThis- linear 
Forchheimer equation was then modified to describe plane 
radial steady gas flow, yielding Equation (44). 

Partial differential equations were developed to 
describe transient gas flow through linear and plane radial 
horizontal systems. 

Numerical solutions were obtained for all four 
boundary conditions considered in the linear system; however, 
time limitations prohibited obtaining similar results for 
the radial system. 

In order to determine whether the solutions Gb= 
tained had converged with respect to the space variable, 
ee of the four sets of boundary conditions specified in 
Table 1 were solved for three different dimensionless space 
increments using that set of dimensionless time increments 
which was considered optimal for the boundary condition in 


question. The three space increments used were 


ieee O20, 0,10, and 0.05, 
The convergence of the solutions with respect to the 
time variable was checked out in-a similar manner. Using 


that dimensionless space increment which was considered 


«i Eo 
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optimal, the four sets of boundary conditions were solved 
employing three different sets of dimensionless time increments. 
The ratio of the dimensionless time increments was maintained 


at two to one, that is 


CA) (At) 2 
DAG a D'B a) iu (54) 
(Ath), (Ath) aL 
where (Ath) as (Ath) pg, and (Ath)o represent the three sets of 


dimensionless time increments that were employed for this 
convergence check. One set, (At,), ris summarized in Table 3. 
From this table it can be seen that the value of the dimen- 
Sionless time increments was increased as the dimensionless 
time increased. 

Each set of the boundary conditions for the linear 
system was solved for dimensionless slip coefficients of 0.0, 
0.2, and 0.4, and for dimensionless inertial resistance co- 
Srfercilents, Oo; U0, 40, 20, 30,740, and 50. The “two sets of 
boundary conditions involving a constant terminal mass pro- 
ducing rate were solved for dimensionless mass flow rates 
(oq,) of -0.10, -0.05, and -0.01. (The negative sign denotes 
production.) For the two boundary conditions involving a 
constant terminal pressure solutions were obtained for. dimen- 
sionless pressures at the producing face of 0.10 and 0.25. 
This resulted in fifty-four solutions each for Cases I and II, 


and thirty-six solutions each for Cases III and IV. 
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Table 3 
Case Summary of the Dimensionless Time Increments (Ath) » 
eee es ee ee ee ee ee ee ee 
li Lc =O 2702 s-7At = 0.04; At = 0.05; 
Bi mg Dy Deas 
At = 0.10; At SU. Oe = 0.30; 
Doan Daou 4 hE, 
Ath = 0.50; Ath = 1.00; Ath = 2.00 
16-20 21-22 23-25 
——— SSNS 
Zeb At = 0.02: At a= 0 504% “At = 0.05; 
D, 33 Dy, fine, 
At =O lTrOr "st = 0.202% AE =" 0.30: 
wa Dig NG 
Ath = 0.50; Ath = 1.00 
16-26 27-35 
eee ee ee 
III Ath = 0.005; Ath = 0.010; Ath = OF 02:3 
1-4 5-8 9-10 
Ath = 0.025; Ath = 0.05; Ath = 0.10; 
11-14 15-18 19-22 
At 20720" At = 0.25% At = 0.50; 
P53 D5 4-27 ges | 
SE = 1.00 
Dyg23g 
————————— 
IV (At, ) = same as for Case III. 
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Graphs for all four boundary conditions have been 
plotted from those results that were obtained at what was 
considered the optimal grid spacing, while these results 
themselves are presented in Appendix E. 

A material balance, as developed in Appendix C, 
was performed in conjunction with these solutions. It 
serves as a check on the numerical procedure used in calcu- 
lating the flow rates at the internal and the external 
boundaries and in calculating the pressure profiles in the 
reservoir. Table 5, which summarizes the material balance 


percentage errors, is presented in the next section. 
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DISCUSSION OF RESULTS 


A. Discussion of Convergence 


Figures 4, 5, 6, and 7 illustrate how variations in 
the dimensionless space increment affected the pressure- 
squared distribution. These curves show that a dimensionless 
Space increment of 0.05 was sufficiently small to ensure con- 
vergence of the solution with respect to the space variable, 

The effect that variations in the dimensionless 
time increment has on the pressure-squared distribution is 
illustrated by Figures 8, 9, 10, and 11. The curves show 
that, except for Case II, the solutions have not quite con- 
verged with respect to the time variable. Since increasing 
the number of time steps increases the computational time, 
it was felt that an optimum between computational accuracy 
and computational time had been reached using the set of 


dimensionless time increments, ( t_) 


D’ AY 


Figures 4 through 11 are all plotted for a dimen- 
sionless inertial resistance coefficient of 30.0 and a di- 
mensionless slip coefficient of 0.2. Since these values are 
the median values of the variables considered, graphs employ- 
ing these values should give a good indication of the conver- 
gence trends for the entire range of variables studied. 

In addition, the following convergence trends were 


observed for the four sets of boundary conditions considered; 
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1. The deviations between solutions obtained using two dif- 
ferent grid spacings decreased as the dimensionless slip 
coefficient increased. 

2. The deviations between solutions obtained using two dif- 
ferent grid spacings increased as the dimensionless iner- 
tial resistance coefficient increased. 

Problems were encountered with the two sets of 
boundary conditions where pressure at the outlet face was 
held constant - for one condition only. When inertial ef- 
fects were negligible (i.e. B = 0.0) but the slippage ef- 
fects were not, the deviations between solutions obtained at 
two different grid spacings were appreciably higher than 
normal. The terms C(u,us) and D(u,us) in Equation (43) be- 
come highly nonlinear when the dimensionless inertial coef- 
ficient is zero. Because slippage effects are not negligible, 
the first-order derivative must also be retained and con- 
vergence problems, which are otherwise negligible, become 


Significant. 
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Figure 4. 
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B. Discussion of Results Obtained at the Optimal Grid Spacing 


Since there is a large number of results included in 
this section, they have been grouped according to boundary 
condition, thus simplifying their labelling. 

Figures 12 to 34 inclusive were obtained employing 
dimensionless space and time increments of 0.05 and (Ath) ar 
respectively. Both increments were considered optimal. The 
primary purpose of these figures is to illustrate the behavior 
of the reservoir with respect to time, and to indicate the 
manner in which the dimensionless inertial resistance and slip 
coefficients affect flow. 

The pressure-squared distributions in the reser- 
voirs, when transient Darcy flow occurs, are indicated by 
Figures 12, 20, 26, and 31. It may be seen that Case I and 
IV stabilized at dimensionless times of 3.0 and 1.0 respec- 
tively, as is indicated by Figures 12 and 31. Cases II and 
III, which depleted at dimensionless times of approximately 
7.5 and 10.0 respectively, are represented by Figures 12 and 
JL. 

Figures 13, 21, and 27, and Table 4 indicate effects 
of -the dimensionless slip coefficient on the pressure- 
squared distributions. From Figure 13 it is obvious that the 
molecular streaming effects become more significant as di- 
mensionless time increases. This is expected since pressure 


in the reservoir decreases with time and the slippage phenome- 
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non becomes more pronounced at lower flowing pressures. It 

is interesting to note that Case I stabilized at approximately 
the same dimensionless time regardless of the slip coefficient; 
however, the stabilized pressure-squared distributions them- 
selves depended considerably upon the value of the dimension- 
less slip coefficient. 

Figure 21 shows that the pressure distributions 
for different dimensionless slip coefficients cross over. 

This can be attributed to the boundary conditions which state 
that the reservoir is being depleted at a constant rate. The 
mean pressure in the reservoir must therefore remain constant 
for any dimensionless time regardless of the Slippage effects. 

Figure 27 indicates the contribution that molecular 
streaming makes on the pressure profiles for Case III. A 
Similar plot for Case IV was not included because the range 
of dimensionless slip coefficients considered did not cause 
any appreciable differences. However, the results obtained 
are summarized in Table 4, 

The effect of the dimensionless inertial resistance 
coefficient on the pressure-squared distribution is exhibited 
by Figures 14, 15, 22, 23, 28, and 32. Deviations from the 
Darcy-continuity solution increased with increasing inertial 
effects. These deviations are more pronounced than the devia- 
tions caused by slippage. When inertial effects are included, 
another resistance to flow is introduced and the pressure- 


squared gradient increases. This steepens the profiles while 
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slippage decreases the resistance to flow and flattens out 
the profiles. Therefore, when both effects are accounted 
for, the corrections for each tend to cancel out. 

Figures 14 and 15 show that the additional pressure 
drop due to the inertial resistance term is more severely 
affected by a decreasing flow rate than the pressure drop 
that is due to the viscous resistance term. The same feature 
can be observed for Case II from Figures 22 and 23. This is 
expected because the pressure drop due to the inertial term 
is dependent on the square of the dimensionless mass flow 
rate while the pressure drop due to the viscous term is 
linearly dependent on the dimensionless mass flow rate. 
Figures 28 and 32, which also show the effect of the dimen- 
Sionless inertial resistance coefficient on the pressure- 
Squared distribution, are self-explanatory. 

The effect of the dimensionless slip coefficient 
and the dimensionless inertial resistance coefficient on the 
pressure at the producing face is illustrated by Figures 16, 
Ppns2Aiprand 8251; These plots are applicable to the constant 
terminal rate cases only - Case I and Case II. It is evident 
from these plots that, for the range of variables considered, 
the inertial effects are again dominant. 

By comparing Figure 16 with Figure 17 and Figure 
24 with Figure 25, it may be observed that a marked difference 
occurs when the mass flow rate is decreased. For a dimension- 


less mass flow rate of -0.01, the curves for different 
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dimensionless inertial resistance coefficients were practically 
superimposed. No conclusions other than the one already 
stated could be inferred from the additional results obtained 
for a flow rate of -0.01. Consequently these results are not 
included. 

Figures 18, 19, 29, 30, 33, and 34 show the effect 
of the dimensionless inertial resistance and slip coefficients 
on the dimensionless mean pressure of the reservoir. 


M 
these figures also represent the depletion history of the 


Since P, = (1.0 - Wo) as indicated in Appendix C, 


specific reservoir in question. 

It is obvious from Figures 18 and 19 that the de- 
pletion rates for Case I increase with increasing inertial 
effects and decrease with increasing molecular streaming. 
This merely verifies the fact that the pressure drop corres- 
ponding to a fixed flow rate increases as the inertial ef- 
fects increase.° In Figure 18 the curves for’B = 30, 40, and 
50 do not extend over the entire range of dimensionless times 
considered since the pressure at the producing face dropped 
below a value of zero. Production was therefore terminated at 
this point. This problem did not occur for a dimensionless 
mass flow rate of -0.05, as evidenced from Figure 19. 

For Case II the mass production rate from the 
reservoir remains constant with respect to time and the de- 
pletion rate remains constant as well. Furthermore, the 


mean pressure is independent of the values of the inertial 
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and the slippage terms because, in Spite of these added ef- 
fects, the net flow from the reservoir is governed by the 
boundary conditions alone. Consequently depletion for all 
values of the inertial and slippage terms followed the same 
straight line decline, and plots similar to Figures 18 and 
19 were not included. 

Figures 29 and 30 and Figures 33 and 34 show that 
the depletion rates for Case III and Case IV decrease with 
increasing inertial effects and increase with increasing 
Slip. This is opposite to the trend in Case I, but never- 
theless, is reasonable for the case of a constant pressure 
at the producing face. This is due to the fact that the 
pressure gradient must increase when the inertial term in- 
creases or when the slippage term decreases. Therefore, 
Since the pressure is fixed at the producing face, the in- 
terior pressures must all be higher. This corresponds to a 
higher mean pressure in the reservoir and thus a lower de- 
pletion rate. The lines corresponding to a dimensionless 
Slip coefficient of 0.4 were not included in Figures 29 and 
30 because they would fuse the family of curves and detract 
from the overall clarity of presentation. For Case IV the 
small differences in the depletion rates due to molecular 
streaming did not warrant including these differences on 


Figures 33 and 34. 
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CASE I 
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Figure 13. 
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Figure 14. Effect of B on Pressure-Squared 


Distribution b = 0.0 (oq, ) a=) 710 





Rigure, L5). 


B 
B 
B 


0.0 —oO——_o— 
20.0—e—— ee — 


40.0—o— —_0— 





Effect of B on Pressure-Squared 


Distribution b 


G80) tego) = -0405 


£20": 
12 ee 7 











a 

pM 0 =O ot - a 
Ooo ce, 2) ee aa ee ——---—-- ee aS ay, 

O 








Jax 





Pagure, 16.) (Bitect of Bi and) b on Outlet Pressure 





(og,) = 0.10 








Outlet Pressure 


es eel 
Os - al i ae yrs a ohh i y 
ook fle va « 
a 
= 


yy i Se "= 
7 ae — 





M 


Ay 


Seige o_o 
Se ee 


b 
b 
b 


0.4 —o— —o— 





18. Effect of B and b on Mean Pressure 


(oq) gene Oe LO 





far 





Figure 19. Effect of B and b on Mean Pressure 
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Figure 22. Effect of B on Pressure-Squared 
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Figure 24.. Effect of B and b on Outlet Pressure 
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Table 4 


B = 10.0 

BP = 0;10 

W 

x 
0.20 0.40 0.60 0.80 
ee eee eee ee ee ee ee ee 
0.0 0.637 0.866 0.947 0.980 

0.632 0.863 0.944 0.979 
0.629 0.861 0.943 0.978 
0.475 0.736 0.870 0.946 
0.470 0,732 0.867 0.944 
0.467 0.729 0.865 0.943 
0.337 0.583 0.757 0.888 
0.331 O.577 0.754 0.886 
0.329 0.575 6 2752 0.885 
On2rak 0.426 0.622 0.812 
O24 Ona le? 0.615 0.809 
0.212 0.415 0.613 0.807 
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C. Discussion of Material Balance Results 
ee ee Pe eDCe RESULTS: 


Table 5 summarizes the errors obtained in the material 
balance calculations which were performed on each of the four 
boundary conditions. The errors, which are of reasonable mag- 
nitude, are markedly smaller for the boundary conditions where 
the flow rates at the boundaries are known and need not be 
calculated. This would suggest that the flow rate calculations 
are the prime source of error in the numerical procedure em- 
ployed. 

Material balance errors will also be introduced by 
approximating the definite integrals by discrete equations, 
as can be seen in Appendix C. 

Problems were once again encountered with Cases III 
and IV, when the inertial effects were set to zero and the 
slippage effects were not. Material balance errors increased 
considerably. This can be attributed to the fact that the 
solutions had not converged, as previously mentioned on page 


42. 
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CONCLUSIONS 


In view of the studies undertaken in this investi- 


gation and the results that were obtained, the following 


conclusions can be formulated: 


18 


The linear Forchheimer equation describing steady hori- 
zontal gas flow through porous media can be derived from 
the Navier-Stokes equation of motion by a procedure simi- 
lar to that employed by Hubbert, but extended to include 
inertial effects. 

The equation thus obtained can be modified to describe 
plane radial flow. 

A mathematical model which accounts for inertial and 
molecular streaming effects can be developed to describe 
transient isothermal gas flow in linear horizontal and 
plane radial porous systems, 

The mathematical model for the finite linear system can 
be solved numerically to yield stable and convergent 
solutions for ihe four sets of boundary conditions that 
were considered. The numerical method employs the Crank- 
Nicholson Implicit procedure, together with a function 
iteration scheme. 

Inertial effects contribute another resistance to gas 
flow through a porous system and therefore increase the 
slopes of the pressure-squared distributions in the 


reservoir. 
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Molecular streaming decreases the flow resistance and 
therefore decreases the slopes of the pressure-squared 
distributions. Consequently the effects tend to cancel 
One another. 

In the study of transient isothermal gas flow through 
porous media, neither the inertial nor the molecular 
streaming effects can be neglected without some former 
investigations to establish that they are indeed negli- 


gible. 
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RECOMMENDATIONS 


The following recommendations are made: 
The results for the finite linear system should be veri- 
fied empirically by running laboratory tests for each of 
the four boundary conditions considered. Care should be 


taken to use homogeneous samples. 


The model describing transient isothermal gas flow through 


a plane radial system should be solved by a numerical 
technique similar to the one employed on the linear hori- 
zontal system. If stability problems are encountered, a 
Backward Difference Implicit procedure, which is less 
accurate but more stable, should be attempted. 

The mathematical model for the two systems should be ex- 
tended to include an absolute permeability which varies 
with distance. This would, however, present problems in 
combining the modified Forchheimer equation with the 
equation of continuity because absolute permeability 
would have to be treated as a variable, 

The procedure could also be extended to include gases 
whose viscosity and compressibility factor vary with 
pressure. A procedure similar to that of Al-Hussainy 


1 (49) 


eta could be employed to allow for these additional 


variables. 



















= s | an ?- 
ormoniicces 
. pala? Si» oe 2) ested " 
rebsm 9%s enotsebnomqost paiwoliot eft 6h. i 
-tzev od blyode medaye taonil odimtt ‘efit x02 etuacront t 
to doss 10% atees viostsxodsl pninags yd ylisoia kc ) bSiz 
ed blyode 9182 .betebLenod anolsibnoo yisbavod wok ed - 
dpuord? wolt esp Ismieritost tnsbensxd patdiizseeb tebom oat 
Isoitemun 56 vd bevioe od biuode meteye {stbsx onaig 5 
~t1iorn xsenil sit no beyolqme sao ez of xsl imie eupindoes 
s ,bsistnuoons eis ameldorq yoilidssa 31 - MBFOES tsaa0s 
aeol et doidw ,siubescqug sata on eoneze8t id | | 
.betqmesis ed bivore veldage stom tud oseuoo8 
~x9 od bluode ameteye owt oft 320d on isoktemedsam of 
esiisvy doidw ytilidsemieq stufoeds 18 shuloai ot bebaes 
nk ameldoxrg jnsesig , revewor binow elaT veoanseib dthy 
edd dtiw sottsype temieddor0® bali thom edd pn! mo: . 
va ilidsomieg stuloeds semsoed ysivaitaos pret: 

-sidsiuey 8 a8 Bedsext ed o2 sad 
seesp shylont o¢ bebnetxe sd oats Biuoo « 
dtiw yisv 1z0f¢pst viilidteasxqmes te pin 
yaiseevh-LA to sted oF sl Lake P 
Isnoltibbs seedt 10% wolis oft 


PiZ 





a 


Both the linear and the 
for infinite reservoirs 
This would exhaust most 
that are encountered in 
studied. 


More accurate numerical 


89 = 


radial systems should be solved 
as well as finite reservoirs. 
of the common boundary conditions 


the two porous systems that were 


integration schemes should be 


employed on the definite integrals in Appendix C in order 


to reduce computational 


errors in the material balance. 
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NOMENCLATURE 
A total cross-sectional area (ne) 
B dimensionless inertial resistance coefficient 
C¥YD,E dimensionless nonlinear terms in PDE 
ES viscous resistance coefficient (oes 
Fy inertial resistance coefficient (L7/) 
Fy inertial forces 
Fy pressure forces 
FL viscous forces 
H total height of sand filter (L) 
Jy hydraulic gradient (F/L) 
K proportionality constant 
L overall length of linear system (L) 
M molecular weight (m) 
L™ ] matrix of nonlinear terms 
eg 8s Shape factors (unitless) 
P pressure (F/L7) 
R universal gas constant (FL/T) 
S surface area ore 
ah temperature (T) 
V volume of the porous system ad 
W mass (m) 
Wy cumulative mass produced (m) 
a vector of known terms in the finite 
difference equation 
b Klinkenberg b-factor (F/L*) 
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b dimensionless slip coefficient 

Cy M/Z eRT (M/FL) 

Cy 1Sa1 C31 Cy specified constants (units vary) 

Cy constant terminal mass flow rate (M/tL*) 

a. characteristic length parameter of porous 
media (L) 

a, mean particle diameter (L) 

g gravitational acceleration (L/t7) 

So conversion factor (L_/Ft*) 

h formation thickness 

hy ,h, fluid heights in manometers terminating 
below and above sand filter 

k absolute (liquid) permeability Ee) 

kK. apparent permeability (17) 

q superficial fluid velocity based on the 
total cross-sectional area (L/t) 

q(t) time-dependent part of (q) 

a radius (L) 

ro radius of a capillary (L) 

t time (t) 

ty dimensionless time 

u dimensionless pressure-squared 

Vv microscopic flow velocity (L/t) 

Ww mass flow rate (M/t) 

Zz gas compressibility factor 

(oq) mass flow rate per total cross-sectional 


area (M/tL@) 
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X,V,Z 


ave 


Subscripts 


X,VrZ 


ave 


Superscripts 


_—_— 


a ae 
distances along the respective rectangular 
coordinates 


proportionality constant that results from 
integration of viscous terms (L7¢) 


proportionality constant. that results from 
integration of inertial terms (L~+) 


mean free path length of a gas (L) 
viscosity (m/Lt) 

kinematic viscosity (L2/t) 
dimensionless length parameter (= ln r) 
density (m/L>) 


porosity (unitless) 


denote the direction in rectangular 
coordinates 


denotes initial conditions 
grid point index along length 
grid point index along time 
denotes standard conditions 
denotes an average value 


denotes reservoir (formation) conditions 


denotes vector quantity 


denotes a dimensionless variable as per 
Equations (B-9) and (B-23) 


iteration number index 
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' Mathematical Notation 


grad..P 
div v 
Dv/Dt 


u—,u—,u 
aex! Trea be|e 


Uz, UT, U 


del operator 

gradient of the scalar (P) 

divergence of the vector (v) 

substantive derivative of the vector (v) 


denote the second-order partial derivatives 
OL. (UU). with respect to x, .G,..& 


denote the first-order partial derivatives 
of (u) with. respect to x, r,.£ 
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Mera PLEIN Dt xX A 


DERIVATION OF THE FORCHHEIMER EQUATIONS 
AND 


THE GENERAL CONTINUITY EQUATION 





A-2 


I. Derivation of the Linear Forchheimer Equation 


The microscopic Navier-Stokes equation of motion in 


the x-direction can be written as: 
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The assumptions inherent in this equation are 
a. Laminar flow 
b. Isothermal flow 
c. Newtonian fluid 
d. The continuous flow channel is saturated with a pure fluid. 
If it is further assumed that 
e. Flow is horizontal 
f. Flow is steady state 
g.-. An applicable equation of state is op = CoP 


then Equation (A-1) can be rewritten as 
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The grouping of the forces permits Equation (A-2a) to be 
written as 


(dF) = (dF_) tt (dF) (A-2b) 


In the following treatment viscosity is considered 
constant and flow is considered to occur through any porous 
media which is homogeneous and isotropic with respect to 
porosity and permeability. It is obvious that Equation (A-2b) 
rewritten in terms of macroscopic variables should consist of 
a pressure-squared gradient which is dependent on two terms, 
namely: 

1. A viscous term which is linearly proportional to a macro- 
scopic flow velocity 

Z. An inertial term which is proportional to the square of 
some macroscopic flow velocity. 

It is possible to define a macroscopic volume ele- 
ment, AV = AxAyAz, which is large enough so that the pro- 
perties of the porous medium do not vary from macroscopic 
element to macroscopic element and small enough so that the 
flow of a fluid through the medium is seen as a continuous 
phenomenon in space. 

A macroscopic fluid velocity in the x-direction 


will be defined according to the following equation: 
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7 * ave (over dAyAz) 


Furthermore, consider” macroscopic flow'to occur parallel to 
the x-axis only. 

The procedure employed to convert the microscopic 
variables involves integrating the terms in the microscopic 
Navier-Stokes equation over the macroscopic volume element 
in order to obtain volumetric averages for the forces in 
question. This method is similar to that employed by M. King 


Hubbert and extends the treatment to include inertial forces. 


a. Pressure Forces 
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(4P,) ye = 7 IO (A-6) 


It should be noted that the pressure-squared gradient 
in Equation (A-6) is the macroscopic pressure-squared gradient 
and is equivalent to the volumetric average of the gradients 
at every point of the fluid contained in the macroscopic volume 


element. 


b... Inertial Terms 
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Consider density to be a function of x only. 
Furthermore, rewrite the derivatives in the above equation 


to include ee Equation (A-7) becomes: 
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The first term in the above equation indicates the 


gradient, in the x-direction, of the microscopic velocity, 


oe and is dependent on three phenomena, namely: 
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1. Expansion and contraction of the flow channels 

2. The change in direction of the flow channels 

3. Pressure drop in the x-direction which causes the fluid 
to expand and the velocity to increase. 

The remaining two terms would be zero if flow was restricted 

to a capillary oriented in the x-direction, because Yee and a 

would then be zero everywhere. However, in a porous medium 

this is the case only when flow is microscopically in the x- 

direction. Everywhere else these two velocities exist. 

Since (av,*/ay) and (av,°/az) have finite values everywhere 

in the macroscopic volume element, the last two terms in 


Equation (A-8) cannot be neglected. The three inertial 


terms will now be treated. 
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A= 7 


This assumption states that the average value of the velocity- 
Squared gradient at every point in the cross-sectional area 
normal to the flow direction is approximately equal to the 
gradient of the average velocity-squared at this cross- 
section. 

In view of Equation (A-4), the following equation 


May be written: 
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The term (oq, / 2x) is the gradient of the macroscopic fluid vel- 
eocity, qe in the direction of flow, and will vary from cross- 
section to cross-section. This term may be approximated by the 


following expression: 


aX. 


where By! and q, are both functions of x. The proportionality 
constant at a certain cross=section, Bi’ accounts for the 
Operator d/ax and must obviously contain dimensions of Cony 
Because the gradient of q,. is negative for production and 
positive for injection, the sign associated with By" will be 


negative and positive respectively. Similarly, because q,, is 


negative for production and positive for injection, the 
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product of the two terms, B,' a, will remain positive for 
both production and injection. It is customary to specify 
the inertial resistance coefficient, which results directly 
from Bi's as positive fom both production and injection. 
Therefore, if By! is specified as positive, an absolute value 
of the macroscopic flow rate must be used in order to retain 
a product which is always positive. As a result, it is pos- 


sible to rewrite Equation (A-13) as: 
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us be B lig | (A-14) 
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OX 
where Bi! is positive for production and injection. 
Employing Equations (A-11), (A-12), and (A-14), 


Equation (A-10) can be rewritten as: 
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Since the macroscopic mass flow rate per total 
cross-sectional area, Pq is constant with respect to x, 


it is possible to rewrite Equation (A-15) as: 
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Integrating with respect to x, the following result 


is obtained: 
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The proportionality constant aay is the average value of By! 
with respect to x, and therefore characterizes the entire 
macroscopic volume element over which integration was per- 


formed. 


Terms 2 and 3: 
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The complexity of the terms © ay bean and 7 = Se 


makes it impossible to rigorously obtain the macroscopic flow 
velocity, qe from these terms even though it is obvious that 
the terms should yield some function of ae In view of this, 


it may be postulated that the following two assumptions are 


reasonable: 
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By virtue of the argument which validated the in- 
clusion of the absolute sign in Equation (A-14), the following 


equations may be written as: 
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Making use of Equations (A-20) and (A-21), Equation 


(A-19) can be written in this form: 
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Because (oq,) is a constant with respect to x, this equation 


can be simplified to read: 
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Combining the three inertial terms into one yields Equation 
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Lumping all three proportionality constants into 
one, Equation (A-25) will become: 
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(dpe. oa (og,) lea, | (A-26) 


The extreme complexity of the proportionality con- 
stant, 8, should be appreciated. It includes the dimensions 
of (ups but the individual proportionality constants that 
comprise it have a length dimension each in a different co- 
OLoinate, direction. (i.e. x, Vv, and.z).. It cannot be written 
in terms of one length parameter, but as some combination of 


all three which will yield dimensions of eae 
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Examining the last two terms first, and considering 
density a function of x only, the following equation 


may be written: 
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As in Equation (A-11) it may be assumed that the 


following statements are valid: 


ae 3 “x ) 
( ) ave (over oAyAz) os sie oven ee) (A-29a) 
dx9d dxdy 

Z d 

hey oacy} 
i= Sp Fenaee __2_avefover $AyAZ) = (n,_oo9p) 
OXdZ dx9Z 


But, because it has been specified that macroscopic 
flow is parallel to the x-axis only, the average values of the 
transversal velocities, vy and vor are zero for any cross- 
section normally to flow. Therefore, it can be concluded 


that (a°v, /axay) and (a°v, /axaz) 


ave(over ¢AyAz) ave(over ¢gAyAz) 


are approximately equal to zero and, as a result, the last 


two terms in Equation (A-27) can be neglected. This agrees 


(41) (42) 


with results obtained by Hubbert and Irmay 


Now consider the remaining three viscous terms. 
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As before, it can be written that: 
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< x ave(over ¢AyAz) 
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3302 ave 3x2 
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> ox 
2 2 
: ox ; (V.-) ave (Over ¢AyAz) 
(eon) a 
5 2 ‘ave 5 2 
Y (wrt oAyAz) Y 
i aq. 
= — 5 (A-32b) 
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2 2 
ic vx) _f ° ve Sve (over ¢AyAzZ) 
32° ave read 922 
(wrt oAyAz) 
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= — 5 (A-32c) 
d dZ 


Using Equation (A-13) as a basis, the following 


approximations may be written: 
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The proportionality constants Oy's to's and O' at 


the particular cross-section being considered account for the 


eo 147 ayli, and S727, respectively, and must 


2 (43) 


operators 3° /3x 
therefore have dimensions of (L “). As stated by Hubbert 
large variations in the velocity, V1 occur over a relatively 
small pore space. The values of On | and a3! are thus ex- 
pected to be large. Variations in cs in the x-direction result 
primarily from the expansion of the fluid due to the pressure 
drop and the value of a5! is expected to be much smaller. 


Therefore, it may be stated that: 

















1 
5 | yong ary (A-34a) 
19° Ix ne! 

' ' 
eBIIK = a aap t (A-34b) 
oa dy a 
4 av. : : 
Consequently, the term 3 owe in Equation (A-30) may be 


neglected. 
Employing Equations (A-32) and (A-33), Equation 


(A-31) may be rewritten as: 


ae Px (2 + a3!) ax 
(dF) ave epg Se Soe a RR Seen Ser gemma (A-35) 
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AX 
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a u Pq.. 
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If the proportionality constants, a, and a2,are 


2 
combined, the average viscous forces for the macroscopic 


volume element may be written as: 


Qa 


(QF) aye 7 TH (eay) (A-36) 
b 


This lumped proportionality constant contains two individual 
proportionality constants which have bey dimensions normal 

to the flow direction (i.e. y and z directions). It suggests 

a reciprocal of the cross-sectional area normal to the macro- 
scopic flow direction. This would definitely imply permeability. 

Employing Equations (A-6), (A-26), and (A-36), the 

microscopic Navier-Stokes equation may be rewritten in terms 

of macroscopic variables which correspond to volumetric aver- 


ages. The result is Equation (A-37) stated below: 


Co 3P- -a B 
poruué—meliias wien ifedshat <7 (oa,) loa, | (A-37) 
> 


2 ox d 


The constants on the right hand side of Equation (A-37) may 


be lumped together in the following way: 


-a 1 

ae oe (A-38a) 
> Kk, 

B 

ri = Fh (A-38b) 


The result will then yield the Forchheimer equation, which 


is written as: 
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9p- 2 u 
-— = — |— (pq.) + F, (oq.) |oq._ | (A-39) 
yx c. k. x b x x 


It is identical to Equation (37). 

The Forchheimer equation thus derived is valid only 
for that macroscopic volume element over which integration 
was performed. However, because the macroscopic volume ele- 
ment was chosen in such a way that it represented the entire 
homogeneous and isotropic porous medium, Equation (A-39) is 
valid for the entire porous system. 

II. Modifying the Forchheimer Equation to Include Plane 
Radial Flow 

It has been established that Equation (A-39) 
describes steady linear gas flow through homogeneous and iso- 
tropic porous media, when inertial effects are accounted for. 
In order to modify this equation to describe plane radial 
flow, coordinate transformations may be used. For horizontal 


flow Equation (A-39) may be written in vector notation as: 


2 u > > > 
-VP“ = — |— (pq) + Fy (og) | oq | (A-40) 
k 
a 


Cc 
Oo 


In rectangular coordinates, the vector quantities become 


, matt vest Se 
VP = lb —t yam + kK — (A-41a) 
ox dy 4 
> > > > 
Gree ee + J dy + k qo (A-41b) 
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For flow in the x-direction only, these expressions reduce to 


2 > ape 
VP = i— (A-42a) 
9x 
gq = i dy. (A-42b) 


In cylindrical coordinates, the same vector quanti- 


ties may be expanded to yield: 


2 s yp Wh yp yp? 
vp’ = § —+ 5 (- sted pa & peer (A-43a) 
ox r 080 OZ 
> > > > 
q = o ae S54 q, + 64, (A-43b) 


If flow is plane radial, the above expressions re- 


duce to: 
5 x ap? 
vp” = § —— (A-44a) 
r 
er 
> > 
Gq = %y 4 (A-44b) 


Since the directions of the x-axis and the r-axis correspond 
for the two flow systems, i and é are identical vectors. 
Therefore, for plane radial flow it is possible to replace 
(x) by (r) everywhere in Equation (A-39). The resulting 


equation that describes such flow is 
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dP 2 U 
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Derivation of the Macroscopic Equation of Continuity 


The derivation of the continuity equation is in- 


Boel’. 


cluded for the sake of completeness and can be accomplished 
by choosing a macroscopic volume element, 


AV, and defining a 
flow vector, Od, as: 


macroscopic mass flow rate 
(eq) 7 = 


(A-46) 
total cross-sectional area 


Summing the flux normal to the surface, over the 


entire surface of this macroscopic element, the following 
equation results ‘442 ; 


—_ > 
Rela eee e= aid .feq)-n ds (A-47) 
AS 
Surface 
where AS is the total surface area. 

By virtue of the Divergence Theorem, 
[f @qy-n as = Jfffv-(pq) av (A-48) 
AS AV 


Consequently the flux leaving the element is equal to the 


negative rate of change of mass within the element, since 
the flux is actually a mass flow rate. 


Therefore: 


f[fa-(eq) av 
AV 


(A-49) 
dt 
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Since (AV) er uag CHAN) Bearer Equation (A-49) 
can be rewritten as 

om 

[{[fv-@qav = - — (A-50) 
AV, /¢ ot 
7 AV. om 

V: (pq) (—— = - — (A-51) 
b a 


AV, is independent of time and may be taken into the deriva- 


Hive, vViz.: 


a in 
ee ie o ——- (——) 
at AV, 
or 
Ve(oq) = -¢ — (A-52) 
at 


which is the general equation of continuity in terms of 


macroscopic variables. 


For linear flow in the horizontal x-direction, 


Equation (A-52) becomes 


9 (pq,) dp 
ge a a =I) — (A-53) 
ox at 


and for plane radial flow Equation (A-52) reduces to 


19 dp 
Sa r (oq) = -> — (A-54) 
r or or 
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DEVELOPMENT OF THE PARTIAL DIFFERENTIAL 
EQUATIONS WHICH ACCOUNT FOR INERTIAL AND 


MOLECULAR STREAMING EFFECTS 


I. 


Linear System 


The following equations describe transient isothermal 


flow of gases through a horizontal linear porous systems 


is 


Quadratic flow equation: 


SOE ait : 
sep eet xn Peereds aa (B-1) 


Continuity equation: 


3 Cea te cry 2 
— (oq) = ae (B=2) 


Slippage correction term: 


b 


ka =Thks¢] e¢ > (B=3) 
Equation of state: 
0 = (ou) P (B= 4) 
ave. 
# cP 


The viscosity of the gas is considered constant and, 


although the gas compressibility factor need not have a value 


of unity (i.e. ideal gases), it is assumed to be constant at 


some average value for the entire flow duration. 


If Equation (B-1) is multiplied by density and com- 


bined with Equation (B-4), a more workable Forchheimer 


equation results: 
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(B-5) 


If Equations (B-3) and (B-5) are combined, the result 


is a quadratic flow equation corrected for molecular stream- 


ing: 


k 


Ope fag 2a IE 
oO 


( sep) (ody) + Fy (04) | oa (B-6) 


In order to extend the range of applicability, Equations 


(B-6) and (B-2) may be rewritten in terms of dimensionless 


parameters. 
aie _ 2 P (oq,) 
is (p+5) 
(pq_) 
P) x ae | 
ox 
where x = x/L 
Bo (B/E. 
b = b/P, 
Lae Wh ea wee 
B = 
rag" 
— 18 
(ody) = a5 
k Pr 
a, = eae 
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These equations may be written as: 


(ody) leas, | 


(B-7) 


(B-8) 


(B-9,a) 
(B-9.b) 


(B-9.c) 


(B-9.d) 


(B-9.e) 


(B-9.f) 
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If common field units are employed, the relationships 


between dimensionless variables and real variables become: 
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where the units in the above equations are: 


u = centipoise 

L = feet 

k = darcies 

Pr = psia 

oe . 3 
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Fy = (feet) 

e = seconds 

«) = fractional 


In order to combine Equations (B-7) and (B-8), the 


7 was adopted, If 


procedure of Tek, Coates, and Katz 
Equation (B-7) is differentiated with respect to the 


variable x, the following equation results: 
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Equation (B-8), a second-order nonlinear partial differen- 
tial equation results: 
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Equation (B=13) can be rewritten more concisely as: 
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Since. the nonlinear terms, C(u,ug) and D(u,usz), contain 
the dimensionless mass flow rate per unit area, Equation 
(B-14) must be solved in conjunction with Equation (B=-7) in 
order to relate (eqy) to the dimensionless pressure 


gradient. 


II. Radial System 


The equations describing transient gas flow through a 


plane radial porous system are: 


1. Quadratic flow equation: 
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3,  Slippage correction term: 
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4, Equation of state: 
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Multiplying Equation (B-15) through by density and elimina- 
ting density from the left hand side, the desired form of 


the Forchheimer equation is obtained, namely: 
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Using the slippage correction term to correct apparent 
permeability for molecular streaming, results in the follow- 


ing expression: 
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Equations (B-16) and (B-20) may be rewritten in terms 


of dimensionless variables as: 
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If the same common field units are employed in the 
radial system as in the linear system, including the well- 
bore radius expressed in feet, the constants relating these 
dimensionless variables to the real variables will also be 
the same as those in Equations (B-10.a), (B-10.b), and 
(B-10.c). 

The same procedure that was employed for the linear 
system may be used to combine Equations (B-21) and (B-22). 
Differentiating Equation (B-21) with respect to the dimen- 


sionless radius yields: 
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Equation (B-28) can be rewritten as: 
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In order to permit closely spaced increments near the 
wellbore and still retain equal increments in the space 


variable, the space variable may be redefined as, €= ln Tos 
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in which case Equation (B-29) may be written as: 
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Since the nonlinear terms C(u,u,), D(u,u,), and 
E(u,u,) contain the flow rate, (oa,), the solution to 


Equation (B-31) must be coupled with the dimensionless 


Forchheimer equation which accounts for molecular stream- 


ing. This equation, in terms of the new independent space 


variable, may be written as: 
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realized that the numerical integration of the describing 


partial differential equation would be more stable if 


Equation (B-22) were not expanded. 


(B-22) should instead be combined into the following 


general form: 
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By employing a finite difference approximation on the space 
variable, a system of nonlinear algebraic equations would 


result, namely: 
ED pO 6 ie ld BN te TY Fa (B-34) 


The properties of the matrix [T] are such that the integrat- 
ing scheme for the above system of equations is more stable 
than the integrating scheme resulting from the discretiza- 
tion of Equation (B=29). In fact if an equation similar 
to Equation (B-33) but in the linear system is discretized, 


the spectral radius of the matrix [T| will be less than 
At 

unity regardless of the value of ( a }. THIS property 
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guarantees stable integration. 
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A material balance was performed by integrating the 
net production from each system with respect to the total 
elapsed time and comparing the dimensionless mean pres- 
sure thus obtained, (Pas with the one obtained by com- 
puting an integrated average of the pressure distribution, 
(Payor at a corresponding time. The difference between 
these two mean pressures was the criterion used to calcu- 


late the material balance percent error, 
ff ) Calculating (Pa 


If it is assumed that the modified gas flow is valid, 
then: 
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a. Linear system 
The initial gas in place may be calculated in the 


following way: 
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b. Radial system 
The initial gas in place is: 
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The mass of gas produced can be calculated in 


the following way: 
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Equations (C-9) and (C-14) are exact expressions for 
the mean pressure that would be obtained in the respective 
porous systems, if the amount of gas produced was used as 
a criterion for the calculation. 

The definite integrals in these two equations were 
evaluated by calculating a two-point average value of the 
dimensionless flow rate for a time increment, and multiply- 
ing this value by the magnitude of the increment. This 
yielded the total mass flow rate out of the reservoir 
during the time increment in question. These flow rates 
were then summed to give a total net production out of the 
reservoir up to some specified time, yielding (Pid 


according to Equation (C-10). 


II. Caiculating (Py) 5 


a. Linear system 


An exact expression for the mean pressure of a linear 
system, as calculated from the pressure profile, is the 


following: 
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Equation (C-15), which was integrated using Simpson's 


Rule, could be approximated by the following discrete 


expression: 
Ax” of = = : 
(Py) 2 Te | Peal + 4 Poe? + 2 Pi =3 ce 
tte oh, err Bare (C-16) 


b. Radial system 


Similarly, it can be written that 
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The definite integral in Equation (C-18) was approxi- 


mated by means of the following discrete equation: 
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Substituting Equation (C-19) into Equation (C-18) results 


in: 
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The percent error in the material balance could then 
be computed for either the linear or the radial system 


according to the following expression: 


% ERROR =.100 __ > (C=-21) 
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DISCRETIZATION OF THE PARTIAL 
DIFFERENTIAL EQUATIONS FOR 


THE LINEAR SYSTEM 





D=-2 


The nonlinear partial differential equation describ- 
ing linear, horizontal, isothermal gas flow through porous 


media is. restated below: 


C(u, us) cers ae D(u, us) sr lis ed (D-1) 


Employing second order correct central-difference approxi- 
mations of the space derivatives, the following discrete 


approximations are obtained: 
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These give rise to an ordinary differential equation, 
represented by Equation (D-4), which approximates the 


partial differential equation: 
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The Crank-Nicholson Implicit procedure uses a forward 
difference approximation on the time derivative, namely: 


du u - u 
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and space derivatives evaluated at a corresponding time 


level of (m+%). Equation (D-4) is therefore rewritten as: 
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Equation (D-6) can be rearranged to yield: 
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This is a general equation for any grid point (%). 
Since an equation of this form must be written for every 
grid point where (U) ntl is unknown, a system of equations 


results, which can be concisely represented by: 


em a Uh te) (D-8) 

The size of this system of nonlinear algebraic equations 
depends upon the number of grid points at which CQ) el must 
be evaluated. As a result the system size depends on the 
set of boundary conditions being considered. Two more 
unknowns are introduced than there are equations. These 
are subsequently removed by use of the two available 
boundary conditions. In order to maintain a consistent 
truncation error, the Neumann boundary conditions were 
approximated by finite difference formulas having the same 
truncation error as the formulas used to discretize the 
differential equation itself (i.e. second-order). 

Each set of boundary conditions yielded a different 
sized matrix with different nonlinear terms. Therefore, 


the boundary conditions are treated separately. 
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CASE I 


14%) At the producing face. 42°=' 1) 


This flow rate is related to the pressure distribution 


through Forchheimer's equation corrected for slippage, 
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Equation (D-7) written for the grid point % = 1 yields the 


following: 
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If Equation (D-12) is used, Equation (D-14) may be 


rewritten as: 
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2. At the external boundary (2 =« LL) 
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Equation (D-7) for the grid point (LL - 1) may be 
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This equation can be rewritten in the following manner: 
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Every other element in the tridiagonal matrix [_M_] and 


the vector (a) are computed using the following general 


equation: 
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The matrix size for Case I will be 


Il = JJ = (LL - 1) (D-20) 


CASE II 


1. At the producing face (2 = 1) 


(og) Tne SF rieksi (0 < m < MM) (D-21) 


Since this boundary condition is identical to the one 
encountered in Case I, the matrix elements M(1,l) and 
M(1,2), and the vector element a(1l) will be defined by 


Equation (D-15). 


2. At the external boundary (yes. LL) 
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When Equation (D-7) is written for the grid point 











2 = LL, the following equation”is»obtained: 
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Making use of the relationship (D-23) and (D-24), 


Equation (D-25) may be written as: 
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At every other grid point, the tridiagonal elements of 
the matrix and the known vector are defined by Equation 


(D-19). 
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The- size of the matrix fer Case II will be: 
Ciy’ = Jd = DG (D-27) 
CASE EILL 
2 At the producing gacefrte = 1) 


Ehatereonderoan (hua si) (D-28) 


Writing Equation (D-7) for the grid point % = 2, the 


following is obtained: 


C2 mts D2, ,m+ 3 2 m+}s 
ig, ae =n (lean es Gad Ve yresseprg 1 oe el 
2 (Ax) 4 (Ax) _ : At, pax | : 





C C 
L L 4, 
: = t mee 43,m+1 ~ a Rs mee asm m| 
2 (Ax) ° 4 (AX) _ ’ 2 (AX) ’ ; ; 
u 
P2 m+} m+ 2,m 
+ e f (D-29) 
4 (AX) 3,m "1 m+ At 
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If the boundary condition dictated by Equation (D-28) 


is introduced Equation (D-29) becomes: 


MUL, oy) o oe Ye a a(l1) (D-30) 


2,m+1 3,m+1 
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1 Cc =< 
where M(1,l1) = : + ae 
At, Cie) 
ie D 
M(1,2)e =nc- eke} . ee FUE! 
2 (Ax) 4 (Ax) | 
Cc 
L 
a(l) = ut fu, onan 
2 (AX) ! 
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2. At the external boundary (Qi=2 LL) 


(Os) tr m+ 


0 (O* < m.<. MM) 


(D-31) 


This boundary condition is identical to the one in Case 


III, and the elements M(II,JJ-1), M(II,JJ), and a(II) are 


governed by Equation (D-26). The remaining elements are 


defined by Equation (D-19), while the matrix size is: 


Idas =apddoxSnat(hh ry ld 


CASE IV 
1. At the producing face (2 = 1) 


= Pp (0 <"m ag MM ) 


(D-32) 


(D-33) 




















a an eed = > i 

7 | aie 

. : bid . 

» er 

Si-a ; 


b A> 27) Yo pee IE ee 
E _ —] = (£,f1)M 9% 
(x6) gt | \ 


(xa) & (xA) 
bt a) 5? 7 ea De, 
(@ g++ uS- u pas “nae 
h as aye (%a) § 
. 
) 7 _ a 
fy S” + Fm $+ sat eee ae 
@? OEE ve hoe s- PI 
= ro sates ¥e each 
sh rashes tenaate 28h-28 
cf 
| re a ” a‘ s 
(f€-a) “(on yer > 0) oO =. “tem, aa! Pa 


? 


9289 ni smo sft of fectsaebi ae > yisbavod 
sis (II)s bas , (UL, II)M panies astbiidte « i a | 
exs ednemels painismes efit .(8S-d) notseuns | 


:ei osie xitstsm sat sdasiiies it oii solamente ts 


(S€-a) (1-1) «= x = rs we _ 
at Pid a area reais ott ¥ 
coe ae EG sa rik ah " 


ea 


a 


(Ls 
thee ane 


(€€~a) | Gm > => ht 
oe ig 

rs hae rat e ‘a 72 “i 

; Na vale ‘ lads Led is) yy J 


" , 
Fearn Vi 


4 





p-r3 


The elements M(1,1), M(1,2) and a(l) are defined by 


Equation (D-30). 


2. At the external boundary te vo) 


Ser dade Ce, CO = MM) (D-34) 
Since this boundary condition is identical to the external 
boundary condition in Case I, the elements M(II,JJ-1), 
M(II,JJ) and a(II) are governed by Equation (D-18). 
Equation (D-19) fixes the remaining elements in the matrix 


equation while the size of the matrix for Case IV is: 
Pet (0. |) ae 92) (D-35) 


Calculating flow rate in the nonlinear terms 


The nonlinear terms C L and D 
2 


cme , contain flow rates, 
a 2 


Q m+ 


(pq,y) Gyr which are related to the dependent variable (u) 
2 


2,m 
through the modified Forchheimer equation. The term (us) 
in this equation was approximated by a finite difference 
approximation with a third-order truncation error. The 
difference equation was written in terms of values at grid 
points which were within the physical system. These finite 
difference equations were obtained by applying the Taylor 


Series expansion (carried out to include the third-order 


derivatives) to the following functions: 
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The derivatives Ceres) and er eer were eliminated by 
selectively combining the eight expansions, the results 
being an expression for (us). 

The following is a summary of the third-order correct 
finite difference equations used to approximate the (us) 


term in the Forchheimer equation: 


a. When 2 L 


27u ~64u, ,,-36u, ,5+288u,,,-215u, 








= 
(uz), = ae - (D-36) 
132 (Ax) 
ie WED s 2 2 See ee LN 2 
-u +6u -3u,.-2u 
+ ~ - 
(uz), = Q+2 Q _ ) Q=-1 (D-37) 
6 (Ax) 
c. When 2= LL-l 
PAE +3u,-6u +u 
+ -1 2- 
(ux) » = ee (D-38) 
6 (Ax) 
Gus When gQ = LL 
215u -288u +36u +64u -27u 
- ~ - -4 
fies) Q g-1 Q-2 g-3 Q (D-39) 


RS 132 (Ax) 
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A copy of the computer program for Case II is in- 
cluded. Since the programs for the other three sets of 
boundary conditions are very similar, copies of these are 
not presented. Therefore, it will suffice to sequentially 
list the differences that do occur, namely: 

1. Different dimensionless time increments. Table 3 dic- 
tates what time increments are to be read into each program. 

2. Different boundary conditions. Statement No. 77 in the 
main line and the subroutine BOUNDS must be altered 
according +0; Table 1. 

3. Different elements in the matrix (M) and the vector (a). 
The subroutine MATRIX must be recalculated, as dictated 
by equations developed in Appendix D. 

4. Different matrix size. Equations in Appendix D state 
the respective matrix sizes. 

5. The production for Cases III and IV need not be termin- 
ated since the pressures in these systems do not drop 
below a value of zero. 

6. The extent: of the loop controlling the flow rate cal- 
culation must be altered so that calculations are 
performed only for those points where flow rates are 


unknown. 
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Computer Program Nomenclature 


Subroutine 


Subroutine 


Subroutine 


Subroutine 


AIN 
AOUT 
AREA 
AC 
3 & 


DC 


BB 


BA 


CAVE 


DAVE 


BOUNDS - sets the boundary conditions. 

MATRIX - calculates the elements of the matrix 
and the known vector. 

SQROOT - inverts the matrix of nonlinear terms 
and multiplies it by the known vector. 
This yields the unknown vector We aye 

MATBAL - calculates the material balance per- 


centage errors at the specified time 


levels. 


matrix of nonlinear terms (M) 

cumulative mass flow into the reservoir 
cumulative mass flow out of the reservoir 
area under the pressure distribution curve 
the terms a, b, c in the quadratic flow 
equation: a(6q,,) * + b(eq,) +c= 0 

the term (b? - 4ac) from the above equation 
dimensionless slip coefficient 
dimensionless inertial resistance coefficient 
vector of known quantities (a) 

nonlinear term C(u,us) 

nonlinear term C(u,us) at time level (m+) 
nonlinear term D(u,u-) 


x 


nonlinear term D(u,us) at time level (m+%) 
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DIFF 


DX 


DTIME 


DT 


DB 


DBB 


El 


EACH 


ERR 


FLOW 
Fl 
F2 
F3 
FS 
Kl 
K2 
KJJ 


IIl=JJI 


LL 


MC 


MB 


E-4 
absolute difference in mean pressure calculated 
by the two ways indicated 
dimensionless space increment 
dimensionless time increment 
dimensionless time increment 
increment in dimensionless slip coefficient 


increment in dimensionless inertial resistance 
coefficient 


error criterion for_convergence of pressure 
distributions (=107>) 


absolute error in pressures for two adjacent 
iterations 


average error in pressures at a grid point for 
two adjacent iterations 


dimensionless flow rate at the producing face 
boundary condition flag (which was not used) 
inertial resistance coefficient flag 

slip coefficient flag 

iteration number flag 

groups of variables used to reduce the length 
of the expressions for the matrix terms 
number of rows and columns in the matrix, 
respectively 

number of grid points at any time level 
number of time steps 


time level at which material balance is to be 
performed 


final time step prior to reservoir depletion 
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Pp 
EP 


PMEAN 
PMEAN1 


PMEAN2 


PCENT 


PAVETI 


PAVEO 


SUM 
TIME 
U 

UU 
UO 
WX 
WWX 


WX1 
WX2 


E=-5 


dimensionless pressure, time level (m+l) 
dimensionless pressure, time level (m) 


dimensionless mean pressure calculated from the 
pressure profiles 


dimensionless mean pressure calculated using 
flow rates 


percent error in material balance 


average flow rate into reservoir for one time 
increment 


average flow rate out of reservoir for one time 
increment 


sum of the individual errors (i.e. }) EACH) 

dimensionless time 

dimensionless pressure-squared at (m+1l) 

dimensionless pressure-squared at (m) 

equivalent to u 
o,m 

(oq,) at time level (m+1) 


(oq,) at time level (m) 


solutions to the quadratic equation 
a (og) 2 + b(oq,) Pie= 0 


equivalent to (u) 


The above summary of nomenclature does not in- 


clude subroutine SQROOT. If this subroutine is treated as 


a black box whose purpose has already been stated, only the 


following variables must be known: 


A 


BA 


the matrix of nonlinear terms (M) 
m+; 

the known vector (a) 

the size of the matrix (M) ths 


pi > 
the solution vector (U) tL 
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REAL UyPyWXoCAVEsDAVEyPPyUUy 1BAsA,C(24,40)yEly ‘EACH, SUM 


Ant rte at 


1,0124 940) »CC (24) 580124) »AC( 24) DC (24) ,ERR(24),2124), 
Z2Bb6eBSeFLUWsTIMEs OTIMEL140) »sPMEAN(40) »pWXLyWX29WWXe9DXaDT, 
3DB8,0D8B 

INTEGER LLyMMyF25F3_,F5,LLLy»MMMygITLI,MC(10) »MB 

COMMON U124940) 9P124540) 9 WX(24 940) »CAVE(24940) sDAVE(24 —_ 
1940) gPP124940) ,UU( 24940) »BA124) ,A124924)4WWX1 24,40) ,DX 
29DT,FLUW 

CUMMON Fisil,MM 




















100 FORMAT(7X,5t10.2) 
1Cl FORMAT(5X,217) 

102 FORMAT(LXs7F10.3) _ 
105 FOCRMAT(/X,516) 

200 FORMAT(* VALUES OF DX,DT,DB.DBB,EL,»,LL,»MM') 
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205 FURMAT(ILH 413Xy"SilP COEFFIB)=",F7.3,* INERTIAL COEFF! 


1,°(88)=",F7.3) 


207 FORMAT(1LHCs35X,*PRESSURE SQUARED DISTRIBUTION") 
210 FORMATILM »45X92FLOWRATE DISTRIBUTION) 
213 FORMAT(30X,"PROBLEM WILL NOT CONVERGE") 

220 FORMATI1L4X,"FIME*,32Xy*DISTANCE(X)= ',35Xy"PMEAN* ) 


AN ANANTH A A AANA RAR RMN Are NA 


221 FORMAT(i4X,'(TD) O2.00 . 0210 0.20 Oe 30 0.40 ! 


i,* 54.50 0.60 C270 0.80 0.90 Ter {PM) *) 


222 FORMATI1L2X,FO.3y,1X%s1lF 7.3) 
223 FORMATI(1X,' oI 


tr OANA LAARARANAA A Oh NANTON Ae ht tc Ae retreat AAs AMAA Q AANA b Bt APRN AA AA AA Atte NA AA ALARA AAA AAA RAH ALARA PLURAL SAA eA NR 


224 FORMATI1L2X_FO.391X,12F 723) 
230 FUORMATA!#—*) 


j _231 FORMAT('1") 


403 FORMAT(*OWXdLsM4#1) TS MULTIVALUED AT M= *,15) 


READ{5,100)0X,DT,DB,0BB,E1 
READS sLOLILL MM I ate iephectieetdnettchaminttinmnitincinenindsipitiibion 


~~ READ(5,105)4MC(1),1=195) 
WRITE(6,200) 


WRITE (6,1 COIDX, OT, DB, UBB EL 


Ce OO” 


WRITE 6,1 01)ILL»MM 
WRITE(6,105)(MC(1),1=1,5) 
READ(5,102)(D07T1ME(M),M=1,MM) 
WRITE(6,102)4D0TIME(M) »M=1,MM) 
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SETTING DAIMENSIGNLESS FLOWRATE 
77 FLOW=—0. 10 


SETTING SLIP COEFFICIENT(B) 


F3=0 

B=-0.19999E+00 
1 F3=F3+1 

B=6+DB 


SETTING INERTIAL COEFFICIENT (38) 


F2=0 
BB=-0.99999E+01 


2 F2=F2tl 
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BB=BB+DBB 
CALL BOUNDS (8,BB) 
LLL=ai~-1 
MMM=MM—1 
TIME=0.0 
DG 99 M=1,MMM 
DT=OTIME(™) 
TIME=TIME+#DTIMEIM) 
C 
C CALCULATING NUNLINEAR COEFFICIENTS (C&D) AT TIMESTEP (M) 
DO14)L=1, LU 
ClisM)=P(LyM)¥(P(LMI+B)/(P (LYM) +4. *BB* (PLL M)4+B) * 
LABS(WX{LyM))) 
Dib yMI=B¥WX(LyM)/ (1 PUL) M) 4B) ¥( PILL Md 44. *BB¥(P(LYM)+B)* 
LABS{WX(LyM)))) a, 
C 
C STORING FIRST GUESS AT (P) AND {WX) AT TIMESTEP {(M+1) 
. PP(L»M+1)=PiiyM) 
| UU(LaM+1)=U(L,M) 
| WWXOL M+ LI=WX(L aM) 
4 CONTINUE 
: | PP(LL+#1,M+1)=PP{LL-1L,M+1) 
UU(LL#1,M+1)=UU(LL-1, M41) 
: wWX (1 »M+] )=FLOW 
¢€ F5=0 
24 CGNTINUE 
Cc 
C CALCULATING THE NONLINEAR COLFFICIENTS (C&D) AT (Mt1) 
DO\5 .L=1,LL 
| Ci »M4+1)=PP1L,M+1)*(PPIL,M41)4+B) /(PPLL, M41 )4+4.*BB*(PP 
14L5M+1)4+8)¥ABS (WWX(LyM#1))) 
Di L yM#L)=B¥wWwX( Ly M4l)/( (PPL L»M+1)4+B) ¥(PP(L, Mtl) +4. *BB* 
LOPPAL »M41)4+B) *ABS(WwX(LyM+1)))) 
CAVELL gy MI=(C1L 4 M)4C(LyM4+1))/2.0 
5 CAVE(LyM)=40(L,M)4+D(L»M4+1))/2.0 
C 
C CCMPUTING SQUARL MATRIX (A) AT TIMESTEP (M41) 
CALL MATRIX (8,BB,™) 
CALL SeROOT {Z) 
Cc 
C SETTING PRESSURES AT (M41) FOR L=1l,Litl 
DO 6 L=lyil 
Uti» MtLd=Z(b) 
IF(UL1,M4+1).G6T.0.005)G0 TO 6 
. . 
€ TERMINATE PRODUCTICN 
DO 69 L=LyLi 
WX(L »M+1)=0.10E-05 
UtLyM+1)=0.005 
69 CONTINUE 
PMEAN(M#1)=0.005 
GO TO 94 
6 P(bLyM4#1)=SQRT(UIL»M4+1)) 
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© CALCULATING THE FLOWRATE AT TIMESTEP (M+1) 
DO 7 L=2,LLL 
IFilsGEsLLLIGO TO 20. ton. 
CCUL)=(-U(L4+29M41)46.¥*U(L41,M+1 )—3. *U(L p M41 )—2.¥*U(L-1, 
1M+1))/46.*DX) 
GO TO 19 

20 CCUL)=12.*U(L+1)Mt1)43. #U(L MF 1)—6.*U(L—1,M41)4UTL—2, 

1M+1))/(6.*DX) 

19 CONTINUE 
BC(L)=-2.0*P(LyM+t1L)/(P(LyM+1)48) 
AC(L)=-4.0*B8 
DCALI=SQRT(BC(L)*BC(L)—-4. O*ACIL) ¥CCIL)) 
IF(CC(L).LE.0.10E-03)G60 TO 17 
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C 
C ALTERNATIVE FOR QUADRATIC CONSTANT = 0.00. 
IF(ABS(AC(L)).GT.0.10E-03)G0 TO 15 
WX(LoM#1)=-1.0%(-CC(LI/BC(L)) 
6G 10:7 
15 WxXl=—(-BC(LI-DC(L))/(2.0*AC(L)) 
WX2=—(-BC(LI40C(L))/(2.0*AC(L)) 
__ TFA (WX1+A8S(WX1)).LT.0510E-04)60 TO 8 posbeade 
LF((WX2+ABS(WX2)).LT.0.10E-04)GU TO 9 
8 WX{LeM+1)=WX1 
LF ( (WX2+ABS(WX2)) »LT.U. LOE-04) WRITE(654039M 
GO'TO\7.° 
17 WX(L»M+i)=0.10E-05 
GOLTSt esi) Se cadericclncchisihivsicicmrnintlowmatnacon Tee 


9 WXiLyMtl)=WX2 
7? CONTINUE 




















. 
C CRECKING TQ See I NONLINEAR TERM HAS CONVERGED 
SUM=0.0 
DC 10 L=ly,LL pe re . 
ERR(LJS=ABS(PP(L »M4+1)—-P(L,M4+1)) 
10 SUM=SUMFERRIL) 
EACH=SUM/FLUAT(ILL) 
DO 13 L=l,iLl 
PP{LyM+1)=P(L,M+1) 
UU{L»M+1)=U(L,M4+1) 
13 WWX{(LsM4+1l)=wxXli,M+1) 
IF{EACH.LE.EL)GO TO li 
F5S=F5+1 
TFIFS.GT.10JGO TO 12 
GO TO 24 
11 CONTINUE 
Cc 
C PROBLEM HAS CONVERGED FOR TIMESTEP (M41) 
C CALCULATING P{MEAN) FROM PRESSURE DISTIBUTION(=PMEANL1) 
PMEAN(M4+1)=0.00 
CG 33 L=2,lLLlL-2 
AREA=0.050* (P(L—-1,M+1)4+4.0*P(iy»M+1)4+P{1L4+1,M41))/3.00 
PMEANIM+1)=PMEANIM4+1L)+AREA 
33 CONTINUE 
GO TQ 99 
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C PROBLEM WItl NCT CONVERGE 
12 wRITE(6,213) 
GO TO 98 
$9 CONTINUE 
¢ 
C wRITING OUT FOR ALL TIMESTEPS 
$4 CONTINUE 
MB=M—1] 
WRITE(6,231) 
WRITE(6,230) 
WRITE( 64230) 
WRITE (6,205)B,8B 
WRITE(6,207) 
WRITE 6,220) es .. 
WRITE{6,221) 
WRITE(6,223) 
TIME=0.00 
DG 78 M=1,MB 
TIME=TIMEF+DTIMETM) 
WRITE16, 224) TIME ,(UIL »M+1) pL=lyliy2) »pPMEAN(M41) 
78 CONTINUE 
WRITE{6,231) 
wRITE(6,230) 
WRITE16,230) 
WRITE(6,205) 0,86 
WRITE(6,210) mw oh 
WRITE(6,223) 
TIME=0.00 
DO 79 M=1,Mb 
TIME=TIME4+DTIME(M) 
WRITE(6, 222) TIME, (WAIL M41) ,L=1,lis2) 
79 CONTINUE 
CALL MATBAL {(MC,DBDTIME,PMEAN) 
C 
C CRECKING EXTENT OCF INERTIAL TERM 
IF(F2.GE.6)GU TO 98 
GG TG 2 
$8 CONTINUE 


C CRECKING EXTENT OF SLIP TERM 
__iF(F3.GE.3)G0 TO 97 
GO TQ 1 
$7 CONTINUE 
STOP 
END 
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SECOND SET UF BUUNDARY CONDITIONS nee 


SUBROUTINE BOUNDS (8,88) 
REAL ‘Us Py yWXy FLOW 


~~ COMMON U(24940) 9P1245940) 9WX(24940) , »CAVE( 24940) »DAVE(24 | 


1940) 9PP(24940) ,UU(24 540) »BA(24)9A1(24924) pWWX1 24940) »DX 
Z29DT,FLOW 


COMMON Fi,it»M™ 


Steet nei rup erate ater vaeteteeeerenn ot rete Saaanony ST DUCC CE UEE STL T SUOTSEDOD 


DG 10 L=lyLL 
P(Lyi)=1.0 
WX{Ly1)=0.106-05 
UtLy1l)=PiLyl)*P(L, 1) 
P(LL+l,1)=P(LL-1,1) 
U{LL#1,1)=U(Li-151). 
DO 13 M=2,M™ 
WX{LL»M)=0.10E-05 
WX1,M)=FLOW 

RETURN 

END 
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SUBRCUTINE MATRIX (8,BB,M) 
REAL UsPyWXeCAVEs DAVE ,yBA,A,U01 24) gK1140) gK2140),KINN140) 5 
LDXeDisByBByPPyUUyWwXx VE 
INTEGER L[isgJJeilsJJJ,FL,iI!I 
COMMGN U(24740) P1124 940) -WX(24540) »CAVE( 24940) »DAVE( 24 
1940) 9 PP( 24540) 9 UU1 24 940) » BAI24) 4 A1 24524) py WWX1 245940) DX 
29DT,FLOW 
COMMON Fl,Liy»M 
LI=LL 
JJ=ELL 
DG 6 J=l,JJ 
| OC 6 i=l,il 
6 Al(i,JI=0.0 
t 
C SETTING FIkST AND LAST ROW OF MATRIX (A) AND KNOWN VECTOR (BA) 
UBIM)=2-0*¥DX¥(2.*P11yM) *WXI1,M)/(10P11)M)4+8)44.*BB*WX(1, 
LM) *ABSTwX01,M)))4002,M) 
K1I{MI=CAVE M)*({U(2,M)—-2.*U(1,.M)+UCIM) {2.*DX¥*¥DX)+ 
LDAVE( 12M) *(U125)M)—-U0IM) )/1(4.*DX)4U(1,M)/0T 
K2(M)=—2.¥*OX¥*¥(CAVE( 19M) /(2.¥*DX*DX)—-DAVE(19M)/14.*DX) )* 
1(12.*PP{ 1 yMtlL)¥WwX1 1 Mtl) /(PP( 1, M41 )4+B)+4.*BBeWX( 19 M4+1)* 
ZABS(WXt(1,M+1))) 
Adis lJ=1.-C0/DTACAVET1L ss MIS(CX*DX) 
Tales DED} 
BACLI=KLIM)—K2{M) 
ALII,JJ—-LI=— CAVE(LL, MIS (0X*DX) 
Alli, JJ J=11.0/DF +CAVEILL »M) /{(DX¥*DX) ) 
KJJIMI=CAVE( Li M)*1(2.0*U{LL—-1,M)—-2.0*U{ LLM) )/1(2.*DX* 
L1DX)+U(Li,MI/DT 
BACLT I)J=KJIIM) 
C 
C SETTING THE REST GF THE DIAGONALS TO APPROPRIATE VALUES 
Ili=II-i oe 
DO 12 I=2,1II 
DO 13 J=laJJ 
IFC ABSTI-J).NE.1)GO TO 13 
IF(J.GT~1)GG TO 7 
AdiyJI=—-CAVEL 1 9M)/12.*DX*UXItDAVE( 1 »MI/(142*DX) 
GG TO 13 
G Ail,JJ=1L-O/UTHLAVEL 1» M)/(DX¥*UX) 
__ GG 10/13 
7 Alig JI=—-CAVE(1,MI/12.0*0DX*DX)-DAVETI»M)/1(4.*DX) 
13 CONTINUE 
BAL ITI=(CAVE(1M)/(2.¥*DX¥*DX) P¥(UCT+1,M)I—2.*U(1»MI+UTI-1 
L1eM))4+(0AVE01 9M) /14.¥*0X))*(U01T4+1,M)-U1I-1,M))4U0T,MI/0T 
lZ2 CONTINUE 
RETURN 
END 
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SUBROUTINE SQROOT (Y) 

REAL A,yBA,AL1 245924) AM(1 24,24) »,AU( 245,24) ,AV124924) 4 
LAA(24 924) 485 (24) 9 V(24) »DIFI24) »PCENTA 24) pXIS 9 VALUE 9 ZX 
29ZXX eo ZV_LZWelQyZRyBO0B 

INTEGER MyKK ye KPeLINeg JP,» JLNGTH, Ky KLy JK oJ 

COMMGN U( 245940) 9P (24540) »pwWX124 540) » CAVE (24940) »DAVEA 24 
1,40) 9PP124940) »UU124,40) »,BA(24) ,A1 24924) yWWX1 24540) DX 
29DT,FLOW 

COMMON Fl,ll,MM a Mare LIER tel mh 

151 FORMAT{LH »29H THE METHOD IS NOT APPLICABLE) 

M=LL 

DO 1 Il=l,™ 

2 AALI,JI=AL1,/) 
1 BB(T)I=BA(1) © 
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C 

v CALCULATICN UF L AND U 
DO 105 K=2 5M 
KK=K-1 


DO 104 J=1,KK 
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AUIK,J)=0.0 
104 AViK,J)=0.0 
105 CONTINUE 
TFL(A(1,1).2EQ20.) GO TO 150 
1C6 X1S=ABS{(A(1_1)) 
AL(1,1)=SQRTIXIS) 
AULislI=Al1,1)/AL(1,1) 
DCG 107 J=2,™M 
AUT 1L,JI=ALLyJISAL(T 1,1) 
107 AL( Je lI=HAldy1LISAUI1y1) 
OC 115 K=2,™ 
KK=K-1 
VALUE=0.0 
CDG 108 J=ly,KK 
108 VALUE=VALUEFALI( Ky J) *AUTJsK) 
ZX=A(KyKI-VALUE 
IF{ZX.EQ.-0.) GO TG 150 
1C9 ZXX=ABS1ZX) 
AL(KsK)I=SGRT(ZXX) 
AU(K,KI=ZX/AL(K,K) 
KP=K+1 
KK=K—-1 
TF(KP.GT.M ) GO TO 115 
DC li2 1=KP,M 
ZV=0.0 
ZW=0.0 
DO 110 LP=1l,KK 
ZV=ZV+AL{KyLP)*AUILP,1) 
110 ZwWw=ZWtAL{1,LP)*AUILP,K) 
AU(Ke1I=(AlKs1I-ZVISALIK,K) 
112 AL(ITsKI=(ALL,»KI-ZWISAULK,K) 
115 CONTINUE 
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C L ANU U ARE CALCULATED 
c PROCEEDING TO CALCULATE L-INVERSE AND U-INVERSE L ee 
DO 119 K=1,M 
119 AV{KyK)I=1.,0/ AUL KK) 
| DO i25 K=2_,™ 
KK=K-1 
DO 122 J=1,KK 
ZQ=0.0 
DO 120 L=JyKK 
120 ZQ=ZQtAL{KyL)*AM(L GJ) 
122 AM{K,J)=-Le/ALTK,K) 
125 CONTINUE 
TJN=M+1 
DO 135 KL=1l,IJN 
K=1LJN-KL 
IF(K.LE.1) GC TO 136 
i2Zé DO 130 JK=1,« 
J=K-JUK 
1FiJ.iT.1) GO TO 135 
iz? ZR=0.0 
JP=J+l1 
_ DO 128 L=JPyK 
128 ZR=ZR#AULJ,L)I*AVIL,K) 
130 AVidsKIH—LZR/AUL I,J) 
135 CONTINUE 
C 
C PROCEEDING TO CALCULATE G-INVEKSE=U-INVERSE*L-INVERSE 
136 DO 140 K=1)™ 
DO 140 J=1,™ 
A{Jr+K)=0.0 
DG 137 L=1l),™ 
137 Ald sKI=AL Se KIFAVI Sy L)FAM( Le K) 
140 CONTINUE 
GO TO 152 
150 wWRITE(65,151) 
152 CONTINUE 
BG 153 1i=1_yM 
Y(ij=0. 
DO 153 J=1l,™ 
153 Yd(IJ=Y(1)+A(1,5)*6AlJ5) 
DO 9 J=1l,M 
BOB=0. 
DO 12 i=ly”™ 
lz BOB=BOB+tAA(J,1)*Y(1) 
DIF(J)=B8(J)-b60B 
PCENT(J)=DIF(J)*100./8B(J) 
9 CONTINUE 
202 RETURN 
ENO 
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SUBRUUTINE MATBAL (MC,yDTIME,PMEAN) 

REAL P,U,WX»)TDyTIME(6) »PMEAN1», PMEAN2, AREA 
- i © 2ePMEANL40) 

INTEGER MyMC(10)9LLyLLLyKol1yMG 


9AIN,AQUT, OUT 
—tatNy REGCUT,REINs DIFF »PCENTIO) »DTIME(40) sPAVEO,PAVET 


COMMON Ul24, 40) 9P 124) 40) » WX 1245940) gCAVE({ 24540) ,VAVE( 24 
1240) »PP(24,40) ,UU( 24,40) ,BA124) 9A( 24524) pWWX1 24540) DX 


220T,FLOW 
COMMON Fl git yMM 


jeotewetonennen mrrrevertetenes 


200 FORMAT(1HOy,26X,*MATERIAL BALANCE CALCULATIONS-% ERROR*) 


201 FORMAT (1HO»21X,"(TD)= 0.100 0.500 1 
| LS 6.000") 
202 FORMAT (24Xy5F8.3) 
11=5 
LS Sel Sc A ae 


C 
C BALANCE ON 5 DIFFERENT DIMENSIONLESS TIMES 


e060 32000', 


DO 20 I=ly,il 


M=MC(1) 


Cc 
C CALCULATING P(MEAN) FROM PRESSURE DISTIBUTION(=PMEANL) 


PMEANL=PMEANIM) 


c 

C CALCULATING PIMEAN) FROM FLOWRATES {(=PMEAWN2) 
AIN=0.00 
AOUT=0.00 
MG=M-1 
DO 15 K=1,MG 
PAVEQ=ABS{(WX11,KItWX(1,K+1))/2.0 
PAVEI=ABS (wX(LL,K)tWXI(LL,K+1))/2.0 
OUT=PAVEO*DTIME(K) 
IN=PAVEL*OTIMELK) 
AGUT=AQUT+OUT 
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AIN=AIN+IN 
15 CONTINUE 
PMEAN2=1.00-ACUT+AIN 
c 
€C CALCULATING PERCENT ERROR IN MATERIAL BALANCE 
60 DIFF=ABS(PMEANI—-PMEAN2) 
PCENT(1)=1LOG.O*DIFF/PMEANL 
20 CONTINUE 
WRITE( 6,200) 
WRITE(6,201) 
WRITE16,202)(PCENT(I),1=1,11) 
RETURN 
END 
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COMPUTER RESULTS 





The computer results that are presented in this 
section were obtained employing the optimal grid spacing, 


namely: 


Ax = 0.05 


At (Ath) a 


Results for all four sets of boundary conditions considered 


are included. 
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CASE I - CONSTANT TERMINAL RATE 
WITH A CONSTANT PRESSURE AT THE 


EXTERNAL BOUNDARY (eq, ).. iG a 
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CASE IV - CONSTANT PRESSURE AT THE 
PRODUCING FACE AND A CONSTANT 
PRESSURE AT THE EXTERNAL BOUNDARY 
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